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1.  Introduction  / 

*  ■  / 

We  have  investigated  a  number  of  topics  related  to  image  / 

restoration,  clutter  reduction,  detection  of  moving  targets,  spatial 
and  temporal  filtering,  and  spectral  estimation.  These  topics 
originated  with  various  aspects  of  the  HALLO  program  but  have  more 
general  applications  including  underwater  detection,  target  identification 
from  undersampled  radar  returns,  deconvolution,  ultrasonics,  detection 
of  points  and  edges,  detection  of  hidden  periodicities  and  speech  processing. 

In  the  course  of  this  investigation,  we  developed  and  applied  the 
following  two  general  methods  of  signal  processing: 

Extrapolation  of  bandlimited  signals  We  introduced  a  numerical 
algorithm  for  estimating  a  signal  f(t)  beyond  a  given  interval  (-T,T) 
under  various  assumptions  involving  mainly  its  spectrum. 

Adaptive  frequency  domain  filtering  and  estimation  We  introduced  a 
method  of  filtering  based  on  the  notion  of  the  running  Fourier  spectrum. 

The  investigation  led  to  a  number  of  reports,  papers  presented  on  various 
conferences  including  a  prize-winning  method  of  spectral  estimation, 
and  papers  published  in  professional  journals.  We  list  below  the  principal 
results  of  the  last  two  years: 

Conference  Papers 

Image  Restoration  (section  2) 

3 

Sixth  Strategic  Space  Symposium  (S  -W),  S.  R.  I.  International, 

Menlo  Park,  California,  March,  1978. 


I 


2 


Adaptive  Extrapolation  and  Hidden  Periodicities  (section  3) 


RADC  Spectral  Estimation  Workshop,  Rome  Air  Development 
Center ,  Rome,  N.  Y,  ,  May  1978.  (The  method  was  applied  to  the 
solution  of  problems  proposed  and  was  a  winner  in  a  spectral  estimation 
competition. 

Improvement  of  Range  Resolution  by  Spectral  Extrapolation 

Third  International  Symposium  on  Ultrasonics,  Imaging,  and  Tissue 
Characterization,  National  Bureau  of  Standards,  Gaither sbure,  Md. 

June  1978.  B 

Adaptive  Frequency  Domain  Estimators 

IEEE  International  Symposium  on  Information  Theory,  Grignano,  Italy, 
June  1979. 


Spectral  Estimation  from  Random  Samples 

IEEE  International  Conference  on  Information  Sciences  and  Systems, 
Patras,  Greece,  July  1979. 

Papers  in  Professional  Journals 

The  Two-to-One  Rule  in  Data  Smoothing 

IEEE  Trans,  on  Information  Theory,  September  1977,  vol.  IT-23 
no. 5,  pp. 631-633. 

Generalized  Sampling  Expansion 

IEEE  Trans,  on  Circuit  &  Systems,  November  1977,  vol.  CAS-24, 
no.  11,  pp.  652-654. 

The  Problem  of  Transmission  Zeros  in  Deconvolution 
IEEE  Trans,  on  Information  Theory,  January  1978,  vol.  IT-24,  no.l. 

pp. 126-128. 

The  Factorization  Problem  for  Time-Limited  F. -notions  and  Trigonometric 
Polynomials 


IEEE  Trans,  on  Circuit  &  Systems,  Jan.  1978,  vol.CAS-25,  no.l. 
pp.  41-45. 


/ 
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Papers  in  Professional  Journals  (cont'd) 

Identification  of  Systems  Driven  by  Non- stationary  Noise 

IEEE  Trans,  on  Information  Theory,  March  1978,  vol.  IT-24,  no.  2, 
pp. 240-244. 

Improvement  of  Range  Resolution  by  Spectral  Extrapolation 

Ultrasonic  Imaging,  vol.  1,  no.  2,  pp.  121-135,  April  1979.  (section  4) 

Detection  of  Hidden  Periodicities  by  Adaptive  Extrapolation 

Acoustics,  Speech  and  Signal  Processing,  vol.  ASSP-27,  no.  5, 
pp.  492-500,  October  1979.  (section  5) 

Theory  and  Applications  of  Running  Transforms 

Submitted  to  the  IEEE  Transactions  on  Circuits  and  Systems. 

Adaptive  Frequency  Domain  Filtering  and  Estimation 

Submitted  to  the  IEEE  Transactions  on  Information  Theory. 

In  the  following  sections,  we  attach  several  of  the  above  papers. 
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Section  2  UNCLASSIFIED 

27.  IMAGE  RESTORATION* 

By:  A.  Papoulis  and  C.  Chamzas+ 


Abstract.  An  extrapolation  method  is  presented  for  recovering 
an  object  from  its  diffraction  limited  image.  The  high  frequency 
components  of  the  object,  lost  due  to  the  finite  size  of  the  aperture, 
are  restored  by  a  numerical  iteration  involving  only  FFT.  The 
method  is  particularly  effective  for  objects  consisting  of  isolated 
points.  Various  illustrations  show  that  the  resulting  resolution  far 
exceeds  the  Rayleigh  limit. 

A  modification  of  the  method  yields  an  algorithm  for  recovering 
a  distant  object  obscured  by  low-frequency  clutter.  The  noise  is 
removed  with  suitable  filtering  and  the  lost  frequency  components  of 
the  object  are  recovered  by  extrapolation. 


1.  Spectral  Extrapolation 

We  shall  develop  a  method  for  restoring  an  object  that  is  distorted 
by  a  diffraction  limited  imaging  system.  For  simplicity,  we  assume 
that  the  object  is  one- dimensional  and  completely  coherent.  The  re¬ 
sults  can  be  readily  extended  to  two- dimensions  and  to  incoherent 
objects. 


In  Fig.  1,  we  show  a  simplified  schematic  of  the  imaging  system. 
We  denote  by  f(x)  and  g0(:c)  the  amplitude  of  the  object  and  of  its  image 
respectively,  and  by  F(u)  and  GQ(u)  their  Fourier  transforms. 
Assuming  that  the  system  is  ideal,  we  have 


<VU> = 


|u|  <r 
W\  >r 


a). 


■where  r  is  the  radius  of  the  aperture. 


♦ - 

Work  sponsored  by  DARPA,  Contract  No.  N00014-67-A- 0438- 0017 

T  A.  Papoulis  is  professor  of  electrical  engineering  at  the  Polytechnic 
Institute  of  New  York,  Farmingdale,  N.  Y.',  11735.  C.  Chamzas  is  a 
Fh.D.  student  at  the  same  school. 
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Object  plane  Fourier  plane  Image  plane 


FIGURE  1  DIFFRACTION- LIMITED  IMAGING  SYSTEM: 

OBJECT  f(x),  IMAGE  gn(x)  AND  FIELD  F(u)  ON 
FOURIER  PLANE  - - .  - 


We  are  given  gQ(x)  and  our  problem  is  to  find  f(x). 
a.  Iteration 

We  shall  show  that,  if  the  object  f(x)  vanishes  for  |x|  >a,  then  it 
can  be  completely  recovered  from  g  (x)  by  numerical  iteration 
(Refs.  1-3).  ° 

For  this  purpose,  we  compute,  first  the  Fourier  transform  G  (u) 
of  gQ(x).  Clearly,  G0(u)  =  F(u)  for  |u|  <r  and  Gq(u)  =  0  for  |u|>r.  0 

W e  truncate  g  (x)  forming  the  function 


fj(x)  = 


gQ(x)  |x|  <a 


We  find  the  Fourier  transform  F^(u)  of  f^(x)  and  form  the  function 

|  Fjfu)  |u|  >r 

Gi(“)=  ,  ,  (3) 

(  F(u)  |u|  <r 


We  find  the  inverse  transform  g,(x)  of  G,(u)  .  This  completes  the 
first  step  of  the  iteration  (Fig.  2). 
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-r  0  r 
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FIGURE  2  FIRST  ITERATION  STARTING  WITH  KNOWN  IMAGE 
.  go(x)  AND  ITS  FOURIER  TRANSFORM  GQ(u) 


The  next  step  is  formed  by  replacing  in  Eq.  (2)  the  function  g  (x) 
by  the  function  g^(x)  and  proceeding  as  in  the  first  step.  ° 

The  nth  step  is  similar:  starting  from  g  .(x),  we  form  the 
function  f  (x)  as  in  Eq.  (2),  the  function  G  (u)as  in  Eq.  (3),  and  its 
inverse  gn(x). 

We  prove  in  Ref.  1  that,  in  the  absence  of  noise,  the  iteration 
converges  to  the  unknown  object  f(x).  This  is  not  the  case  if  noise  is 
present.  However,  even  then  the  method  can  be  used  to  enhance  the 
distorted  image.  This  is  done  by  terminating  the  iteration  at  an 
appropriate  step  so  as  to  minimize  the  overall  error  (Ref.  1). 

b.  Point  objects  .  If  additional  information  about  the  object  is 
available,  then  the  iteration  speed  increases  and  the  noise  problem 
can  be  relaxed.  A  case  of  particular  interest  in  many  applications 
involves  objects  consisting  of  distinct  points  (neighboring  stars,  for 
example).  In  this  case,  f(x)  consists  of  sharp  peaks  at  certain  points 
and  the  problem  is  to  find  the  location  and  the  amplitude  of  these  peaks. 
The  number  of  points  need  not  be  known  in  advance  although  this 
knowledge  improves  the  convergence  of  the  iteration. 
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Adaptive  extrapolation.  To  determine  f(x),  we  start  the  iteration 
as  iiTTTgTZT  Since  t(x)  consists  of  impulses,  the  function  gn(x)  ex¬ 
hibits  maxima  for  sufficiently  large  n.  When  this  is  observed,  a 
threshold  level  is  set  and  only  the  region  of  gn(x)  above  that  level  is 
retained.  The  remaining  part,  below  the  threshold  level,  is 
replaced  by  zero  yielding  the  function  f  (x)  shown  in  Fig.  3.  The 
iteration  step  is  completed  as  before.  Tne  threshold  level  is  set  low 
at  first  and  is  raised  as  the  iteration  progresses.  As  we  show  in  the 
following  illustrations,  the  object  f(x)  can  be  recovered  even  in  the 
presence  of  considerable  noise. 


FIGURE  3  TRUNCATION  OF  g  (x)  BELOW  A  THRESHOLD  LEVEL 
YIELDING  fn+1(x)  n 


Example  1.  In  Fig.  4,  we  show  an  object  f(x)  consisting  of  two  points 
in  a  noisy  background,  and  its  image  gD(x).  It  is  evident  from  the 
figure  that  the  two  points  are  completely  merged.  The  results  of  the 
iteration  are  shown  for  n=4  and  n=30.  As  we  see,  the  amplitudes  and 
locations  of  the  two  points  are  recovered  for  n=30. 

Example  2.  In  Fig.  5,  we  show  an  object  consisting  of  6  points  in  a 
noisy  background  and  its  image  g  (x).  At  the  12th  iteration,  we 
observe  6  maxima  but  the  amplitudes  are  not  correct.  The  6  points 
are  essentially  recovered  at  the  100th  iteration. 
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w  (x),  we  apply  the  preceding  method  suitably  modified:  Reversing 
tl?e  role  of  high  and  low  frequencies,  we  can  show  that  the  iteration, 
starting  from  wq(x)  yields  f(x).  The  details  are  omitted. 


FIGURE  5  f(x):  UNKNOWN  OBJECT  CONSISTING  OF  6  POINTS 

CONTAMINATED  BY  NOISE,  gD(x):  KNOWN  IMAGE, 
(x):  RESULT  OF  12TH  ITERATION, 

100TH  ITERATION 


Example  4.  In  Fig.  7,  we  show  the  unknown  object  f(x)  and  the  date 
y(x)  containing  the  noise  n(x).  Filtering  of  y(x)  yields  the  signal 
wo(x).  As  we  see  from  the  figure,  the  noise  is  eliminated  but  the 
signal  f(x)  is  strongly  distorted.  To  find  f(x),  we  apply  the  iteration 
method.  At  the  10th  step,  we  obtain  the  signal wj~(x)  shown  in  Fig. 7  . 
This  is  essentially  identical  to  the  unknown  object  f(x). 


Sioo(x):  ^SULT 
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FIGURE  6  (a)  UNKNOWN  OBJECT,  (b)  KNOWN  IMAGE 
.(c)  4TH  ITERATION,  (d)  24TH  ITERATION 


FIGURE  7  f  (x);  UNKNOWN  OBJECT,  y  (x):  IMAGE  OF  f(x)  VIEWED 
THROUGH  CLUTTER,  w0(x):  FILTER  OUTPUT, 
w3(x):  3RD  ITERATION,  w1Q(x):  10TH  ITERATION 
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Section  3  . 

ADAPTIVE  EXTRAPOLATION  AND  HIDDEN  PERIODICITIES 


A.  PAPOULIS  and  C.  CHAMZAS 


Department  of  Electrical  Engineering 
Polytechnic  Institute  of  New  York 
Route  110,  Farmingdale,  N.  Y.  11735 


1.  Introduction 


An  important  problem  in  many  applications  is  the  determination  of  the 
frequency  components  of  a  signal 


(1) 


in  terms  of  the  segment 


w^t) 


f(t)  +  n(t) 
0 


It  |  <  T 
It  |  >  T 


(2) 


of  f(t)  containing  the  noise  component  n(t).  The  signal  f(t)  is  not  known  for 
every  t  for  a  variety  of  reasons: 

The  signal  f(t)  can  be  written  as  a  sum  of  exponentials  for  a  limited 
time  only  (Voice;  non-stationary  processes). 

The  available  time  of  observation  is  limited  (sun  spots;  weather  trends). 

Measurements  are  limited  by  instrument  constraints  (Michelson  inter¬ 
ferometer;  band-limited  channels). 

The  unknown  frequencies  a),  and  coefficients  c.  can  be  determined 
simply  with  ordinary  Fourier  transforms  if  the  time  of  observation  2T  is 
large  compared  to  all  the  periods  T.  =  2n/(JUi  .  This  is  not,  however,  the 
case  if  T  is  of  the  order  T^  particularly  if  the  noise  component  of  w^(t)  is 
not  negligible.  In  this  paper,  we  present  a  method  which,  as  we  hope  to 


-13 


show,  is  reliable  even  in  such  extreme  cases. 

The  method  involves  only  FFT  and  it  is  based  on  earlier  results 
dealing  with  the  problem  of  extrapolating  band-limited  functions  Ll,  2] 
We  review  for  easy  reference  the  relevant  parts  of  these  results. 

2.  Extrapolation  of  band-limited  functions 

Consider  a  function  f(t)  with  Fourier  transform  F(ui)  such  that 


F(uu)  =  0  |uu|  >  o 


We  form  the  function 


Wj(t)  = 


It  I  <  T 
It  I  >T 


obtained  by  truncating  f(t)  as  in  Fig.  1.  We  shall  determine  f(t)  in  terms  of 
w^(t)  by  numerical  iteration. 

First  step.  We  compute  the  Fourier  transform  Wjiu)  of  w1(t),  form  the 
function 


FjOu)  = 


iW^(ui)  |oj|<a 

jo  |uu|  > a 


Compute  the  inverse  transform  f^(t)  of  F^(w),  and  form  the  function 


w2(t)  = 


vx(t)  =  f(t)  |t  |  <  T 
,(t)  |t  |  > T 


and  its  Fourier  transform  W^fuJ). 

This  completes  the  first  step  of  the  iteration  (Fig.l). 


2 
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Hence, 

E  -E 
n  n+1 

-  2  2 

-  +  TS  j  |F(»)-Wn+I(,»)|  d» 

jt  I  <  T  |uj|>cr 


(10) 


In  Cl]  and  [2]  we  show  that  f  (t)-f(t)  as  n~».  This  is  not  true  if  the 
given  segment  w^t)  of  f(t)  is  noisy  as  in  (2).  In  this  case,  a  satisfactory 
estimate  of  f(t)  can  be  found  by  early  termination  of  the  iteration.  [2] 

Note.  From  (10)  it  follows  that  the  mean-square  error  E_  is  a  monoton 
decreasing  function  and  since  it  is  positive,  it  tends  to  a  limit.  This  does  ret 
prove  the  convergence  of  (9)  because  the  limit  need  not  be  zero.  It  shows, 
however,  that 

1 

]  En"En+r  °  n~” 

j  Hence, 


J  [f(t)-f 

|  It  |<T 

I 

I 
I 

Althoughthe  functions  f(t)  and  f  n  (t )  are  band-limited,  (11)  does  not  imply 
that  f(t)  -*fn(t)  because  there  is  no  lower  bound  on  the  energy  concentration 
of  band-limited  functions  in  a  finite  interval  [l,  3]  .  For  example,  the 
prolated  spheroidal  functions  cp  (t)  are  band-limited,  their  energy  equals 
one  but  their  energy  concentration  in  the  interval  (-T,  T)  tends  to  zero  as 
n-®.  This  is  the  case  because  the  eigenvalues  X  of  the  underlying 
l  integral  equation  tend  to  zero  as  n~®. 

We  mention  without  elaboration  that,  in  the  discrete  version  of  the 
problem,  the  convergence  of  the  iteration  can  be  deduced  from  (11)  under 
suitable  conditions.  The  reason  is  that  the  corresponding  eigenvalues  are 
finitely  many,  therefore,  they  have  a  positive  minimum. 


»mJ 


dt 


(11) 
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3.  Adaptive  extrapolation 

The  preceding  method  was  based  on  the  assumption  that  the  unknown 
function  f(t)  is  band-limited.  This  information  was  used  to  reduce  the 
error  in  the  estimation  of  f(t)  twice  at  each  iteration  step.  The  speed  of 
iteration  can  be  increased  and  the  effects  of  noise  can  be  reduced  if 
additional  a  priori  information  about  f(t)  is  available.  Suppose,  for 
example,  that  the  size  of  the  band  of  1  (tu)  is  known  but  its  precise  location 
•  j-s  unknown.  We  then  choose  a  constant  a  sufficiently  large  for  F(uu)  to 
vanish  outside  the  integral  (-a,  a)  and  proceed  as  in  Sec.  2.  As  the 
■  iteration  progresses,  the  form  of  W  (uj)  suggests  appropriate  reduction  of 
the  assumed  band  of  f(t). 

The  adaptive  extrapolation  method  is  particularly  effective  if  f(t)  is  a 
sum  of  exponentials  as  in  (1).  In  this  case,  F(uu)  consists  of  impulses 
(lines)  as  in  Fig.  2; 


m 

F(u>)  =  2tt  )  c.dOu-u;.)  (12) 

i=l 


and  our  problem  is  to  determine  their  locations  oi.  and  amplitudes  c.  in 
terms  of  the  known  segment  wn(t)  of  f(t).  1  1 


To  solve  this  problem,  we  select  a  constant 


a  larger  than  the  largest 

•'  W  I, 


possible  value  of  ti>.  and  we  proceed  with  the  iteration  until  W  (uu)  take ^ 
significant  values  only  in  a  subset  Bn  of  the  band  (-a,  a)  of  f(tf  (Fig.  3). 
This  suggests  that  the  unknown  frequencies  are  in  B  .  When  this  is 
observed,  the  function  Th(uj)  of  the  nth  iteration  step^s  obtained  from  the 
following  modification  of  (7): 


i 


F  (IB)  = 
n 


(Fi&3)In  the  above,  B  is  the  set  of  points  such  that  W  (uo)  exceeds  a 
threshold  level  e  .  n  n 


uu  e  B 

n 


uu  e  B 

n 


W  (UU) 
nv  ' 


uue  B 


uu  e  B 


(13) 


(14) 
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and  Bn  its  complement  .  The  process  is  repeated  until  W  (tu)  approach 
the  unknown  spectrum.  This  can  be  checked  by  ccmparingnthe  inverse 
w  (t)  of  Wn(uu)  with  the  known  segment  of  f(t). 

Notes.  As  the  following  examples  show,  the  unknown  frequencies  can  be 
found  even  if  the  data  are  noisy  and  the  constant  T  is  smaller  than  the 
smallest  period  1\  . 

The  choice  of  the  threshold  level  e  is  dictated  by  two  conflicting 
factors:  For  a  speedy  convergence  ancf  reduction  of  the  noise  component, 
e  must  be  large.  It  must  be  sufficiently  small  so  that  all  frequency 
components  of  F(iu)  are  in  the  set  B  .  Thus  e  is  small  at  first  and  it 
increases  as  the  iteration  progresses. 

The  accuracy  of  the  method  depends  on  the  number  m  of  the  unknown 
components  and  their  relative  locations  and  amplitudes.  If  some 
components  are  small  compared  to  the  maximum  c.,  it  is  possible  that  they 
could  be  lost.  However,  if  the  noise  is  sufficiently  small,  they  can  be 
recovered  by  substracting  the  significant  components  and  repeating  the 
iteration. 

A  priori  knowledge  of  the  number  m  of  the  unknown  frequencies  is 
useful  but  not  essential. 

If,  at  the  nth  iteration  step,  all  frequency  components  of  f (t )  are  in  the 
set  Bn  ,  then  the  resulting  mean- square  error  reduction  is  given  by  (10) 
mutatis  mutandis. 


4.  Illustrations 


We  conclude  with  a  digital  implementation  of  the  above  method.  The 
computations  were  performed  with  a  PDP  11  minicomputer  (single  precision) 
and  the  FFT  size  was  N=  256.  The  known  segment  w^(t)  of  f(t)  contains  the 
first  30  points.  In  figure  4  the  unknown  signal  consists  of  3  cosine  waves. 
Determination  of  their  amplitudes  and  frequencies  is  not  apparent  from  the 
Fourier  transform  W^(iu),  figure  4b,  of  the  given  segment.  The  resolution 
is  improved  considerably  within  a  few  steps  of  the  iteration.  In  the  12th 
step  the  three  cosine  terms  have  been  revealed,  figure  4c,  and  in  the  63rd 
step  complete  estimation  of  frequencies  and  amplitudes  have  been  achieved, 
(figure  4d). 

In  the  next  example,  figure  5,  the  same  signal  is  used  but  we  added 
10%  white  noise  uniformly  distributed.  The  signal  is  again  recovered 
completely  in  100  steps.  In  figure  6  the  unknown  signal  consisting  of  two 
exponentials  has  been  corrupted  by  20%  white  noise.  To  recover  the  signal 
the  iteration  needs  30  steps. 
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Thus,  as  we  see  from  these  illustrations  the  unknown  frequencies  can 
be  found  even  if  the  data  are  noisy  and  the  constant  T  is  smaller  than  the 
smallest  period  in  f(t). 
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FIGURE  4.  (a)  F(uj):  The  Fourier  transform  of  the  noiseless  unknown 

signal. 

(b)  W^('jj):  The  Fourier  transform  of  the  given  segment  W^t). 

(c) ,  (d):  The  result  of  the  12th  and  63rd  iteration 

(£n  indicates  fcbe  threshold  level  at  the  nth  step) 
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(a)  F(m).  The  Fourier  transform  of  the  unknown  signal  'f(t). 

(b)  The  Fourier  transform  of  the  given  segment  w^(t) 

(c) ,  (d):  The  result  of  the  12th  and  100th  iteration. 


FIGURE  6.  (a)  F(uu):  Unknown  signal  consisting  of  two  impulses 

contaminated  by  noise. 

(b)  W^(ai):  Fourier  transform  of  the  known  segment  w^(t) 

(c) ,(d):  The  result  of  the  4th  and  30th  iteration. 
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Section  4 

IMPROVEMENT  OF  RANGE  RESOLUTION  BY 
SPECTRAL  EXTRAPOLATION 

A.  Papoulis'  and  C.  Chamzas 

Polytechnic  Institute  of  New  York 
Department  of  Electrical  Engineering 
Farmingdale,  New  York  11735 


Under  various  simplifying  assumptions,  the  reflected 
signal  y(t)  in  the  interrogation  of  a  substance  by  an  ultrasonic 
wave  is  a  convolution  of  the  transmitted  signal  x(t)  with  a 
function  h(t)  that  is  related  to  the  reflection  coefficient  of  the 
medium  in  the  direction  of  propagation.  The  function  h(t)  can, 
in  principle,  be  determined  by  deconvolution.  However,  since 
the  band  B  of  the  spectrum  X(oj)  of  x(t)  is  finite,  the  frequency 
components  of  h(t)  outside  B  cannot  be  found  reliably.  In  this 
paper,  a  method  is  presented  for  extrapolating  X(u)  beyond  B. 

The  resulting  increase  in  resolution  is  limited  only  by  the  level 
of  noise.  The  method  is  particularly  effective  if  'n(t)  is  a  sum 
of  impulses. 

Keywords:  Deconvolution:  diagnostics;  echoes;  extrapolation; 
resolution;  ultrasonic. 

1.  Introduction 

Ultrasonic  waves  are  used  in  metallurgy,  in  medicine,  and  in 
other  areas  to  determine  the  structure  of  a  medium.  The  medium 
is  interrogated  by  a  narrow  beam  (fig.  1)  and  the  reflected  signal 
y(t)  is  used  to  determine  various  properties  of  the  medium,  in 
particular,  the  location  of  its  surface  of  discontinuity.  Under 
various  simplifying  assumptions,  the  signal  y(t)  can  be  expressed 
as  a  convolution  integral: 


yit)  =  /  x(t-r)h(T)dT  (1) 

-m 

where  x(t)  is  the  transmitted  signal  (fig.  2a)  and  h(t)  is  a  function 
related  to  the  reflection  coefficients  of  the  medium  in  the  direc¬ 
tion  of  propagation.  The  variable  t  is  proportional  to  the  distance 
along  the  beam.  The  assumptions  leading  to  eq.  (I)  and  the  rela¬ 
tionship  between  h(t)  and  the  parameters  of  the  medium  will  not  be 
considered  here. 

If  the  medium  consists  of  homogeneous  layers,  then  h(t)  is  a 
sum  of  impulses  as  in  figure  2b: 
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h(t)  3  /  c.  5(t-t.)  (2) 

i 

where  the  points  correspond  to  the  locations  of  the  separation 
surfaces,  and  the  coefficients  ci  are  related  to  the  reflection  co 
efficients  at  these  surfaces.  In  this  case,  eq.  (1)  yields 


y(t)  =  T  cj  x<t-ti)  (3) 

i 

and  the  problem  is  to  determine  the  unknowns  c.  and  t.  in  terms 
of  the  observed  signal  y(t) . 


If  x(t)  is  a  pulse  whose  duration  d  is  smaller  than  the  mini¬ 
mum  distance  tj-tj  between  neighboring  impulses,  then  the  vari¬ 
ous  terms  in  eq.  (3)  do  not  overlap  (fig.  2c),  hence,  the  constants 
ci  and  ti  can  be  readily  found.  This  is  not  the  case,  however, 
if  d  is  larger  than  ti-tj  (fig.  2d).  The  resolution  of  the  system, 
i.  e.  ,  the  smallest  distance  between  reflecting  surfaces  that  can  be 
detected  without  special  processing  is  thus,  proportional  to  the 
duration  d  of  x(t).  To  improve  the  resolution,  we  must  decrease 
d,  or,  equivalently,  we  must  design  a  transducer  with  large  band¬ 
width.  In  this  paper,  we  present  a  method  for  increasing  the  reso¬ 
lution  of  a  given  system  by  numerical  processing  of  the  signal 
y(t). 


In  principle,  h(t)  can  be  determined  simply  in  terms  of  x(t) 
and  y(t)  by  deconvolution.  Indeed,  taking  Fourier  transforms  of 
both  sides  of  eq.  (1),  we  conclude  from  the  convolution  theorem 
that 


This  yields  the  Fourier  transform  H(w)  of  the  unknown  h(t)  in 
terms  of  the  Fourier  transforms  X(u)  and  Y(w)  oi  the  signals  x(t) 
and  y(t)  respectively.  To  find  h(t),  it  suffices,  therefore,  to  com¬ 
pute  the  inverse  transform  of  the  ratio  Y(u»)/X(u). 

However,  the  resulting  values  of  H(cu)  in  the  region  of  the  fre¬ 
quency  axis  where  X(w)  is  of  the  order  of  the  background  noise  are 
not  reliable.  For  a  typical  transducer,  X(u)  is  a  curve  as  in  fig¬ 
ure  3  taking  significant  values  in  the  band  (ui.uZ)  only;  hence, 

H(<«>)  can  be  determined  reliably  only  in  this  band.  Our  problem, 
therefore,  is  to  extrapolate  H(<«>)  for  <*»  <  wl  and  w  >  uZ.  We  shall 
do  so  under  the  assumption  that  h(t)  is  a  sum  of  impulses  as  in 
eq.  (2).  We  note  that,  since  X(u)  is  not  strictly  bandlimited,  the 
frequencies  ul  and  uZ  are  rather  arbitrary.  If  they  are  so  chosen 
that  X(<*>)  is  significantly  larger  than  the  noise  level,  then  the  er¬ 
ror  in  the  estimation  of  H(u)  by  the  ratio  Y(a))/X(u)  is  small.  How¬ 
ever,  the  extrapolation  problem  is  then  more  difficult  because  the 
length  wZ-wl  of  the  resulting  interval  in  which  H(u)  is  assumed 
known  is  small. 

We  note,  finally,  that  in  the  presence  of  noise,  the  best  esti¬ 
mate  of  H(u)  in  the  interval  (wl,«Z)  is  not  the  ratio  Y(h)/X(u).  It 
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Fig.  3  (a)  Typical  transmitted  signal  x(t) 

(b)  Spectrum  of  x(t) . 


can  be  shown  [4]  that  the  estimate  of  H(<j)  in  the  vicinity  of  the  low 
values  of  X(o)  can  be  improved  if  the  properties  of  the  noise  are 
known. 

The  proposed  method  of  determining  h(t)  in  terms  of  only  a 
segment  of  its  Fourier  transform  H(u),  is  a  nonlinear  adaptive 
modification  of  an  extrapolation  method  developed  in  reference 
[2].  In  the  next  section,  we  review  briefly  the  relevant  parts  of 
this  paper. 


2.  Spectral  extrapolation  of  time-limited  signals 

Consider  a  function  h(t)  with  Fourier  transform  H(<u).  We 
assume  that  only  the  segment 


WjM  = 


H(u) 

0 


u  e  B 


u  5  B 


(5) 


of  H(u)  is  known,  where  B  is  a  certain  band  on  the  frequency  axis 
and  Bis  its  complement.  In  general.  h(t)  cannot  be  determined  in 
terms  of  W|(w).  However,  if  h(t)  is  time-limited,  i.  e.  .  if 

h(t)  =  0  for  1 1 1  >  T  (6) 

then  [21  it  can  be  uniquely  determined  in  terms  of  W^ui). 

In  the  ultrasonics  problem,  the  duration  2T  of  h(t)  is  pro¬ 
portional  to  the  thickness  of  the  medium,  and  the  band  B  is  the 
interval  (ui,ti>2)  in  which  X(w)  takes  significant  values.  The  func¬ 
tion  Wj(u)  equals  the  ratio  Y(u)/X(w)  for  every  u  in  B. 
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Using  a  numerical  iteration  involving  only  FFT's,  we  shall 
determine  H(<*»)  for  every  w. 

First  step.  We  compute  the  inverse  transform 


Wjlt)  dw 

B 

of  Wj(u).  We  form  the  function 


ht(t)  = 


Wj(t) 

0 


M  <  T 


|t|  >  T 


(?) 


(8) 


obtained  by  truncating  Wj(t)  for  |t|  >  T  (fig.  4).  We  compute  the 
Fourier  transform 

T 

Hj(u)  =  /  h1(t)e”j“t  dt  (9) 

-T 


of  hj(t),  and  form  the  function 

j  Wj  (u)  =  H(u)  ueB 


W2(«)  =  J 
nth  iteration. 


Hj(u)  u  e  5 

We  compute  the  inverse  transform 


(10) 


wn(t)  = 


"Sr 


f  W  (<i))e^nl‘,t  du> 


(11) 


of  the  function  Wn(u)  determined  at  the  end  of  the  n-1  iteration. 
We  form  the  function 

j  w„(t>  |t|  <  T 

\(t)  =  J  (12) 

I  0  |t|  >T 

obtained  by  truncating  wn(t)  for  |t|  >  T  (fig.  5).  We  compute 
the  Fourier  transform 


T 

Hn(">  =  fh n(t)e’j“t  dt 


(13) 
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Fig.  5  Description  of  nth  iteration 


of  h  (t),  and  form  the  function 
n 

Wn+l(u)  =  PMH  M  +  [1  -  P(u)]Hn(«)  (14) 

where  PM  is  a  rectangular  window: 

1  u  e  3 

PM  =  (15) 

0 

as  in  figure  5. 

We  have  shown  [2,  31  that,  in  the  absence  of  noise,  the 
sequence  hn(t)  so  formed  tends  to  the  unknown  signal  h(t)  as  n 
tends  to  infinity.  The  effects  of  noise  are  considered  in  refer¬ 
ence  [2]  . 

In  eq.  (14),  we  used  for  the  term  P(u)  a  rectangular  window. 
This,  however,  is  not  necessary  [5,  6]  .  It  can  be  shown  that  the 
iteration  converges  for  any  P(u)  provided  that 


0  <  PM  <  1 


(16) 
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We  mention  also  that  the  speed  of  convergence  increases 
if  eq.  (14)  is  replaced  by  the  equation 


HM 


A  H  M 
n  n 


u  «  B 
u  e  B 


(17) 


where  A  is  a  suitable  factor. 

Q 

3.  Adaptive  extrapolation 

We  now  assume  that  h(t)  is  different  from  zero  not  in  the 
entire  interval  (-T.  T)  but  only  in  a  subset  b  of  this  interval. 

Thus,  h(t)  =  0  for  every  t  in  the  complement  b  of  b.  In  the  prob¬ 
lem  under  consideration,  b  is  a  set  consisting  of  a  finite  num  - 
ber  of  points  t}  [see  eq.  (3)]  .  If  the  set  b  is  known,  then  we  de¬ 
termine  the  function  hn(t)  of  the  nth  iteration  such  that 

Iw  (t)  t  s  b 

n  .  (18) 

0  t  s  b 

This  is  a  modified  form  of  eq.  (12)  and  it  results  in  a  reduction 
of  the  effect  of  noise  and  an  increase  in  the  speed  of  conver¬ 
gence. 

In  our  problem,  the  set  b  is  not  known.  In  fact,  our  objec¬ 
tive  is  to  find  it.  However,  the  information  that  it  consists  of 
a  finite  number  of  points  can  be  used  to  reduce  the  size  of  the 
truncation  interval.  This  is  done  as  follows: 

We  start  the  iteration  as  in  Sec.  2.  As  the  iteration  pro¬ 
gresses,  the  function  wn(t)  takes  significant  values  only  in  a 
subset  of  the  interval  (-T.  T)  because  wn(t)  tends  to  a  sum  of 
impulses  as  n  -  •  .  This  set  is  denoted  by  bn  and  is  so  defined 
that 

wn(t)  >  en  t  e  bQ 

wn(t>  <sn  tsin  (19) 

where  en  is  a  suitable  threshold  level.  It  is  desirable  that  en 
be  as  large  as  possible  subject  to  the  condition  that  all  points 
tj  be  included  in  the  corresponding  set  bn.  We  have  found  from 
a  number  of  numerical  calculations  that  a  reasonable  choice  is 
the  minimum 

eQ  =  minflw^jd)  |  ]  1  e  bn-l 

of  the  signal  w^  ^(t)  in  the  set  bn-l  of  the  preceding  iteration 
step. 

We  next  form  the  function  hn(t)  by  truncating  wn(t)  in  the  com¬ 
plement  bn  (fig.  6): 
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Fig.  6  Adaptive  reduction  of  truncation  set;  s  ;  threshold  level. 

n 


ha(t> 


W  (t) 
n 

0 


t  e  b 


t  e  b 


(21) 


we  compute  its  Fourier  transform  Hn(u),  and  continue  as  in 
eq.  (14).  A  nested  sequence  of  seta  b 


b  .  3  b  3  b  .  3 
n-1  —  n  —  n+1  — 


results  which  for  reasonably  small  noise  levels  converges  to  the 
set  of  points  t*  .  If.  however,  some  of  the  coefficients  cj  of  h(t) 
[see  eq.  (2)]  are  small,  the  corresponding  points  tj  might  not 
be  in  all  the  sets  bn  -  This  depends  on  the  level  of  noise  but 
precise  conditions  cannot  be  easily  established  in  general.  We 
are  in  the  process  of  determining  sufficient  conditions  for  the 
complete  recovery  of  h(t)  for  various  special  cases  involving  a 
small  number  of  points  t.. 

The  following  modification  of  the  method  results  in  a  re¬ 
duction  of  the  effects  of  noise  and  round-off  errors:  the  set 
ba  does  not  necessarily  change  at  each  iteration  step.  i.  e.  . 
it  is  possible  that  bn_ l  =  bn.  We  continue  the  iteration  until 
bn  is  a  proper  subset  of  bn_j.  When  this  is  observed,  we  re¬ 
place  the  signal  wn(t)  by  the  signal  wi(t)  of  the  first  iteration 
step,  and  start  again  where  now  the  interval  (-T.  T)  is  replaced 
by  the  set  bn.  This  reduces  the  total  energy  of  the  noise,  and 
it  eliminates  the  accumulation  of  errors  from  the  preceding 
steps. 

4.  Numerical  illustration 


We  shall  use  the  method  to  determine  the  location  and  the 
reflection  coefficients  of  a  synthetic  medium  consisting  of  six 
homogeneous  layers.  The  corresponding  impulse  response  is 
an  impulse  train  consisting  of  six  unequal  impulses  as  in  fig¬ 
ure  7.  The  observed  signal  is  the  sum 

z(t)  =  y(t)  +  n(t)  (22) 

shown  in  figure  8a.  The  component  y(t)  is  obtained  by  con¬ 
volving  h(t)  with  the  signal  x(t)  of  fig.  8b. 
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Fig.  7  Unknown  impulse  respons  1 


iii»  n(ti 


-(2.8  -4.4  -T  0.0  T  6.4  (2.8 
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Fig.  8  (a)  Observed  signal  z(t)  =  y(t)  +  n(t)  consisting  of 

reflected  signal  y(t)  and  noise  component  n(t); 
(b)  Transmitted  signal  x(t) 


6 

y(t)  =  x(t)  *  h(t)  =  £  c.  x(t-t.)  (23) 

i=  1 

The  component  n(t)  is  white  noise  with  uniform  distribution, 
zero  mean  and  standard  derivation  <7=0.  05.  The  signal  x(t)  is 
the  output  of  a  physical  transducer  (see  also  fig.  3),  and  it  de¬ 
termines  the  scale  of  the  time  axis.  All  other  signals  are  com¬ 
puter  generated. 

Our  objective  is  to  determine  numerically  h(t)  in  terms  of 
z(t)  and  x(t).  The  computations  are  carried  out  digitally  on  a 
single-precision  minicomputer  using  as  sampling  interval 
5  =  0.  lys.  It  is  assumed  that  the  unknown  impulses  are  in  the 
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interval  (-T.  T)  where  T  =  4us.  This  interval  is  generally 
known  from  the  physical  description  of  the  problem  (outer  sur¬ 
faces  of  interrogated  medium),  or  it  can  be  determined  from 
the  observed  data.  In  our  case,  the  interval  2T  contains  80 
sample  points.  All  Fourier  transforms  are  computed  with  an 
FFT  size  N  =  256.  The  total  processing  interval  is,  thus, 

N6  =  25.  6us,  that  is,  about  three  times  the  truncation  interval 
2T.  The  interval  N6  must  be  sufficiently  large  to  avoid  alias¬ 
ing  errors.  The  sampling  interval  5  determines  the  limit  of 
resolution.  The  parameters  N  and  6  are  not  critical. 

Since  it  is  known  that  the  segment  of  z(t)  outside  the  inter¬ 
val  (-T,  T)  is  noise,  it  is  truncated  prior  to  any  processing. 

We  shall  first  attempt  to  determine  h(t)  by  direct  decon¬ 
volution.  It  will  become  evident  from  the  computations  that 
the  results  do  not  give  an  adequate  estimate  of  h(t).  For  this 
purpose,  we  compute  the  Fourier  transforms  Z(u)  and  X(iu)  of 
z(t)  and  x(t)  respectively  and  form  the  ratio 


HaM 


ZM 


(24) 


(fig.  9a).  The  inverse  transform  ha(t)  of  Ha(u)  is  the  result 
of  direct  deconvolution  (fig.  9b).  The  unknown  impulses  are  not 
recognizable. 

The  effect  of  the  noise  is  most  pronounced  outside  the  in¬ 
terval  (uj.cjj).  To  reduce  it,  we  truncate  Ha(u>)  outside  this 
interval.  The  resulting  function  W j (u)  and  its  inverse  wj  (t) 
are  shown  in  figure  10.  Again,  wj  (t)  is  not  a  satisfactory  es¬ 
timator  of  h(t)  because  the  frequencies  outside  the  interval 
(<■>},  big)  are  eliminated. 

To  determine  the  missing  frequencies  of  H(w),  we  apply 
the  iteration  described  in  Sec.  3  starting  with  the  function  Wj(w) 
of  figure  10.  The  results  of  the  iteration  are  shown  in  figure 
11a  for  various  values  of  n.  In  figure  lib,  we  show  in  detail 
the  function  wn(t)  for  n  =  40.  In  figure  12  ,  we  plot  the  nor¬ 
malized  mean-square  error 

255  255 

ea  =  |  )  [h(k6)  -  wo(k6)12  E  =  7h2(k6)  (25) 

k=0  k=0 

as  a  function  of  n.  As  we  see  from  these  results,  the  unknown 
h(t)  is  essentially  fully  recovered  for  n  =  40.  In  fact,  the  lo¬ 
cations  tj  of  the  six  impulses  are  located  exactly  (see  figure 


We  next  repeat  the  iteration  using  as  the  window  not  the 
pulse  P(«)  of  eq.  (15),  but  the  function 


P,(<4  =  £  |  X(u)  | 


(26) 
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Fig.  9  Estimate  of  h(t)  obtained  through  direct  deconvolution 

(a)  Ratio  =  Z(u)/X(u)  of  spectra  of  z(t)  and  x(t) 

(b)  Inverse  transform  h_(t)  of  H.(u). 


Fig.  10  Estimate  of  h(t)  obtained  through  deconvolution  and 
truncation  (a)  Reliable  segment  Wj(u)  of 
(b)  Inverse  transform  w^(t)  of  Wj(w) 


where  1c  equals  the  maximum  of  |  X(u)  | .  This  function  is  be¬ 
tween  zero  and  one  and  is  chosen  because  it  favors  the  fre¬ 
quencies  for  which  |X(u)|  is  large,  i.  e.  ,  the  portion  of  the 
frequency  axis  where  Ha(u)  is  reliably  determined  [see  eq. 
(24)]  .  The  window  Pl(o)),  however,  has  the  drawback  that  it 
distorts  the  component  of  Z(u>)  due  to  the  useful  signal  y(t). 

In  figure  1  3a,  we  show  the  results  of  the  iteration.  The 
function  w4o(t)  is  shown  in  figure  13b  and  the  normalized  mean 
square  error  en  in  figure  12  .  As  we  see  from  those  figures, 
the  error  is  reduced  rapidly  at  the  first  15  steps,  but  it  re- 
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(a)  tb) 


Fig.  11  (a)  Results  of  adaptive  iteration  with  rectangular  window 
for  n  =  8,  16,  24,  32,  and  40.  (b)  Detailed  plot  of 
w  (t)  for  n  =  40. 

Q 


Fig.  12  Mean-square  error 


N-l 

*a*H  Ch(k6)  -wa(k6)]2 
krO 

»;  With  rectangular  window,  b:  With  window  Pj(w)  pro¬ 
portional  to  |  X(u)  | 
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Fig.  13  (a)  Results  of  adaptive  iteration  with  window  P^(u>)  for 
ns  8,  16,  24,  32,  40.  (b)  Detailed  plot  for  n  =  40. 


mains  essentially  constant  for  n  >  20.  The  reason  is  that  in 
the  early  stages  of  the  iteration,  the  distortion  due  to  the  noise 
is  rapidly  reduced,  and  what  remains  is  the  distortion  of  the 
useful  component  of  Z(w)  due  to  the  fact  that  P[(<j)  is  not  con¬ 
stant  in  truncation  interval.  It  seems,  therefore,  that  a  com¬ 
bination  of  the  two  windows  leads  to  better  results:  We  start 
with  the  window  ^(<u)  until  the  set  bn  of  the  iteration  is  only  a 
small  part  of  the  interval  (-T,  T).  We  then  continue  with  the 
rectangular  window  of  eq.  (15). 

We  conclude  with  an  observation  concerning  the  meaning 
of  the  term  •resolution'’.  This  term  is  used  extensively  in 
several  areas,  and  it  has  been  given  various  definitions. 

These  definitions,  however,  lead  to  a  measure  of  resolution 
that  is  of  the  order  of  the  duration  d  of  the  transmitted  signal 
x(t),  or,  equivalently,  inversely  proportional  to  the  band¬ 
width  of  x(t).  (In  optics,  the  resulting  number  is  inversely 
proportional  to  the  diameter  of  the  aperture.  )  The  measure 
so  described  is  adequate  if  we  rely  on  the  direct  observation 
of  the  data  for  the  determination  -ji  the  minimum  spacing  be¬ 
tween  impulses  (point  objects,  in  optics).  It  is  not,  however, 
useful  if  the  data  are  subjected  to  elaborate  processing.  As 
we  have  seen,  if  the  data  contain  no  noise,  then  we  can  sepa¬ 
rate  two  or  more  arbitrarily  close  points  no  matter  how  large 
d  is.  The  resolution  is  limited  only  by  the  size  of  the  sam¬ 
pling  interval  6.  This  is  not  so  if  noise  is  present;  however, 
in  this  case  also  the  limitation  is  not  3implv  the  size  of  d.  but 
it  depends  on  the  level  of  the  noise  and  on  the  form  of  the  un¬ 
known  signal  h(t). 
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Section  5 

Detection  of  Hidden  Periodicities  by 
Adaptive  Extrapolation 

ATHANASIOS  PAPOUUS,  fellow.  IEEE,  AND  CHRISTODOULOS  CHAMZAS 


Abstract— A  method  is  presented  for  determining  the  hsimonic  com¬ 
ponents  of  a  noisy  signal  by  nonlinear  extrapolation  beyond  the  data 
interval.  The  method  is  based  on  an  algorithm  that  adaptively  reduces 
(he  spectral  components  due  to  noise. 

I.  Introduction 

AN  important  problem  in  many  applications  is  the  deter¬ 
mination  of  the  frequency  components  of  a  signal 
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/(0  =  £  W 

i*  1 


/«*»ft 


in  terms  of  the  segment  (data) 


W|  (f)  =* 


/(f)  +  n(f) 
0 


ui<r 

ifi>r 


(i) 


(2) 


of  /(f)  containing  the  noise  component  n(f).  The  data  are 
known  for  If  I  <  T  only  for  a  variety  of  reasons: 

1)  The  signal  /(f)  can  be  written  as  a  sum  of  exponentials 
for  a  limited  time  only  (voice;  nonstationary  processes). 

2)  The  available  time  of  observation  is  limited  (sun  spots; 
weather  trends). 

3)  Measurements  are  limited  by  instrument  constraints 
(Michelson  interferometer;  diffraction-limited  imaging). 
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Fig.  1.  (a)  The  unknown  signal  /(/)  and  iu  Fourier  transform  F( ui). 
(b)  First  iteration  starting  with  known  segment  wjU). 


A 


—  Wn  (“t 

—  'n  «-> 


JSL 


-<r  a  ui 


Fig.  1.  nth  iteration. 


The  unknown  frequencies  ut  and  coefficients  c,  can  be 
determined  simply  with  ordinary  Fourier  transforms  if  the 
time  of  observation  IT  is  large  compared  to  all  the  periods- 
r,  *  2ir/to<  and  their  differences.  This  is  not,  however,  the 
case  if  T  is  of  the  order  Tt  -  T,,  particularly  if  the  noise  com¬ 
ponent  n(t)  is  not  negligible.  In  this  paper  we  present  a 
method  which,  as  we  hope  to  show,  is  reliable  even  if  T  is 
small  and  the  data  are  noisy. 

The  method  involves  only  FFT  and  it  is  based  on  earlier 
results  dealing  with  the  problem  of  extrapolating  band-limited 
functions  (1],  [2],  We  review  (for  easy  reference)  the  relevant 
parts  of  these  results. 

II.  Extrapolation  of  Band-Limited  Functions 
Consider  a  function  /(f)  with  the  Fourier  transform  F(cj) 
such  that 


F,(  «)• 


w,M 

o 


M<o 

I  to  I  >  a. 


(5) 


We  compute  the  inverse  transform  /t(f)  of  F.(to),  and  form 
the  function 


w2(f) 


w,(r)  =/(f) 

"./i(0 


ni<r 

u\>r 


(6) 


and  its  Fourier  transform  W2(to). 

This  completes  the  First  Step  of  the  iteration  (Fig.  1). 
nth  Step:  We  form  the  function  (Fig.  2) 


/■„(<*>) 


IV„(cj) 

.0 


|ui|  <  a 
|w|  >  cr 


(7) 


<c(to)*0  |to|>o. 

We  form  the  function 


w,(f) 


J/(o 

lo 


ui<r 

ifi>r 


(3) 

(4) 


obtained  by  truncatuig  /(f)  as  in  Fig.  1.  We  shall  determine 
/(f)  in  terms  of  w,(f)  by  numerical  iteration. 

First  Step:  We  compute  the  Fourier  transform  (V, (cj|  of 
w,(f)  and  form  the  function 


where  (Vn(co)  is  the  function  obtained  at  the  end  of  the  pre¬ 
ceding  step  and  compute  the  inverse  transform  /„(f)  of  Fn(ut). 
We  form  the  function 


i(0  = 


7(0 

7,(0 


ui<r 
ifi  >  r 


(8) 


and  compute  its  Fourier  transform  (to). 

If  /(f)  is  approximated  by  /„(0,  the  resulting  mean-square 
error  is  given  by 
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r  *  1  Ca 

£n  mJ  lf(0-Ur)l2  Jf—J  IF(cj)  -  Fn(u)l*  du. 

(9) 

We  maintain  that  this  error  decreases  twice  at  each  iteration 
step.  Indeed, 


]2df+  / 

*'in>  r- 


[/(0-/„(r))l<it 


I  i m-uoi1 

"'Irl  <  T 

But  [see  (7)  and  (8)| 

f  m  -  /„(')  1 1  dt  «  f  [/(r)  -  w„M(r)]ldf 

-in>r 

=■  7-  I  IF(to)  -  (V„..l(cj)|2  dco 
-it--.. 

I  *  ^♦i(w)!2  icj 

**  •'I  W|  >  <X 

i  r 

I  IF(cj)  -  Jcj. 

And  since  the  last  integral  [see  (9)] ,  we  obtain 

£W„-,  ■  I  [/(0-/.WI1* 

-'ui<  r 

+  7- I  IK'i.-iMI2 


(10) 


because  /•Yw)  =  0  for  |cj|  >  a. 

In  [  1  ]  and  [2]  we  show  that  /„(r)  —/(f)  as  n  —  <®.  This  is 
not  true  if  the  given  segment  w,(f)  of  /(f)  is  noisy  as  in  (2). 
In  this  case,  a  satisfactory  estimate  of  /(f)  can  be  found  by 
early  termination  of  the  iteration  [2j. 

Note:  From  (10)  it  follows  that  the  mean-square  error  En 
is  a  monoton  decreasing  function  and  smce  it  is  positive  it 
tends  to  a  limit.  This  does  not  prove  the  convergence  of  (9) 
because  the  limit  need  not  be  zero.  It  shows,  however,  that 

E„'E„,  i-0 
Hence. 


I  [/(f)-/„(f)|2dr-*0 


-  UK  T 


(ID 


.Although  the  functions /(f)  and /„(f)  are  band  limited,  (1 1) 
does  not  imply  that  /(f)  — /„(f)  because  there  is  no  lower 
bound  on  the  energy  concentration  of  band-limited  functions 
in  a  finite  interval  [1).  [3).  For  example,  the  prolate  spheroi¬ 
dal  functions  p „(f)  are  band  limited;  their  energy  equals  one 
but  their  energy  concentration  in  the  interval  (-T,T)  tends 
to  zero  as  n  —  °®.  This  is  the  case  because  the  eigenvalues 
X,  of  the  underlying  integral  equation  tend  to  zero  as  n  — 

We  mention  without  elaboration  that,  in  the  discrete  version 
of  the  problem,  the  convergence  of  the  iteration  can  be  de¬ 
duced  from  (11)  under  suitable  conditions.  The  reason  is 
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1 

Fig.  3.  Fourier  transform  of  the  unknown  signal. 
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Fig.  4.  Truncation  of  W*tw)  below  a  threshold  level  «„  yielding  Fniw>. 


that  the  corresponding  eigenvalues  are  finitely  many,  there¬ 
fore,  they  'nave  a  positive  minimum  [4] . 


III.  Adaptive  Extrapolation 
The  preceding  method  was  based  on  the  assumption  that 
the  unknown  function  f(t)  is  band  limited.  This  informa¬ 
tion  was  used  to  reduce  the  error  in  the  estimation  of  /(f) 
twice  at  each  iteration  step.  The  speed  of  iteration  can  be 
increased  and  the  effects  of  noise  can  be  reduced  if  addi¬ 
tional  a  priori  information  about  /(f)  is  available.  Suppose, 
for  example,  that  the  size  of  the  band  of  F( oj)  is  known  but 
its  precise  location  is  unknown.  We  then  choose  a  constant 
o  such  that  F( u)  vanishes  outside  the  integral  (-o.  o)  and 
proceed  as  in  Section  II.  As  die  iteration  progresses,  the 
form  of  W„(<o)  suggests  appropriate  reduction  of  the  as¬ 
sumed  band  of  /(f). 

The  adaptive  extrapolation  method  is  particularly  effective 
if  /(f)  is  a  sum  of  exponentials  as  in  (1).  In  this  case,  /"(cj) 
consists  of  impulses  (lines)  as  in  Fig.  3: 


F(u j)  =  2»r  £  C(S(tj  -  o;(), 

i-i 


(12) 


and  our  problem  is  to  determine  their  locations  ui,  and  ampli¬ 
tudes  c ( in  terms  of  the  known  segment  wt(r)  of  /(f). 

To  solve  this  problem,  we  select  a  constant  o  larger  than 
the  largest  possible  value  of  and  we  proceed  with  the 
iteration  until  (V„(cj)  takes  significant  values  only  in  a  sub¬ 
set  B„  of  the  band  (-o,  a)  of  /(f)  (Fig.  4).  This  suggests  that 
the  unknown  frequencies  are  in  Bn.  When  this  is  observed, 
the  function  Fn(io)  of  the  nth  iteration  step  is  obtained  from 
the  following  modification  of  (7) 


En(cj) 


KO*)  uea, 

0  wSfl, 


(13) 


(Fig.  4)  where  B„  is  the  complement  of  Bn.  The  process  is 
repeated  with  occasional  reduction  of  the  size  of  Bn  as  further 
evidence  suggests,  and  it  terminates  when  w„(f)  is  essentially 
a  sum  of  exponentials.  Another  application  of  the  method 
is  discussed  in  [5]  in  the  context  of  deconvolution. 
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Discussion 

The  adaptive  extrapolation  method  is  essentially  empirical. 
Although,  as  we  see  in  the  following  examples,  it  works  well 
in  a  number  of  cases,  there  is  no  a  priori  certainty  that  in  a 
given  problem  it  will  converge  to  the  unknown  signal.  In  fact, 
if  some  of  the  components  c,  of  f(t)  are  relatively  small,  they 
might  be  lost. 

The  accuracy  and  reliability  of  the  method  depends  on  a 
number  of  parameters:  total  number  of  unknown  frequencies, 
possibly  prior  knowledge  of  this  number,  relative  sizes  of  am¬ 
plitudes  c,  and  frequencies  to,,  noise  level,  length  2  T  of  the 
data  interval,  and  available  FFT  size  .V.  A  precise  statement, 
even  empirical,  of  the  importance  of  all  these  factors  cannot 
be  made:  it  would  depend  on  many  parameters.  We  are  in 
the  process  of  determining,  empirically,  the  limits  of  the 
method  for  a  number  of  special  cases.  We  comment  below, 
briefly,  on  certain  empirical  criteria  for  selecting  the  set  B„ 
and  on  the  limitations  due  to  sampling. 

For  the  subset  5,  introduced  in  (13)  we  select  the  set  of 
points  such  that  the  magnitude  of  W,(cj)  exceeds  a  threshold 
level e„: 

lWn(ou)l>e,  u  Sfl„ 

(14) 

|W,(cj)|<e,  u£S,. 

The  choice  of  e„  is  dictated  by  two  conflicting  require¬ 
ments:  for  a  speedy  convergence  and  noise  reduction.  e„  must 
be  large;  it  must  be  sufficiently  small  so  that  all  frequency 
components  of  f(r)  remain  in  B„.  In  the  examples  given  be¬ 
low  we  used  the  following  method  for  determining  e„. 

We  first  find  the  minimum  Mn. ,  of  !Wn.,(co)t  in  the  set 

5„-,: 

Mn.,  -min  (15) 

If  e„-,  is  greater  than  t ,  where  u  is  a  constant  less 

than  one,  then  we  do  not  change  the  threshold  level.  Ife«_, 
is  less  than  uMn  . ,  then  we  choose  e„  *  i xMn  _ , .  Thus, 

en  *max {e„„,,Milfft-i}-  (16) 

In  the  examples,  u  is  chosen  between  0.9  and  0.99. 
jVumencal  Considerations 

The  numerical  implementation  of  the  method  involves  the 
discrete  signals 

F„=F(nu0) 

obtained  by  sampling /(f)  and  F(u> ). 

Suppose,  first,  that  the  problem  is  inherently  discrete,  i.e., 
that  we  wish  to  find  the  spectrum  of  a  sequence  /„  from  in¬ 
complete  data.  Clearly,  the  discrete  version  of  the  iteration 
and  of  the  band-limited  assumption  are  self-evident.  How¬ 
ever.  the  assumption  that  /(f)  is  a  sum  of  sine  waves  has  no 
obvious  discrete  version.  It  corresponds,  loosely,  to  the  as¬ 
sumption  that  the  smallest  distance  of  the  nonzero  frequencies 
is  large  compared  to  one  (no  ''neighboring  frequencies"  are 
present).  If  this  is  the  case,  then  the  unknown  frequencies 
can  be  determined  exactly,  provided  that  the  data  interval  is 
not  too  small  and  the  noise  level  is  reasonable. 
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Fi*.  5.  Discrete  spectrum  IFmI  of  fn  ■  «  0.3,1V"  256). 


We  tum  now  to  our  main  problem:  the  numerical  deter¬ 
mination  of  the  frequencies  of  an  analog  signal.  We  assume 
that  the  FFT  size  iV  is  specified.  It  suffices,  therefore,  to 
select  the  size  ta  of  the  sampling  interval.  As  we  know  [  1 J , 
the  frequency  interval  is  then  determined  because  * 
2nlNta.  Since  the  data  interval  is  2T.  the  number  M  of 
available  samples  equals  2 T/t0.  The  choice  of  M  is  guided 
by  the  following  considerations:  if  M  «.V,  then  the  itera¬ 
tion  might  converge  to  the  wrong  frequencies.  If  M  is  large, 
then  the  abasing  errors  are  large. 

It  appears  from  our  experience  that  ,Vf  » ,V/4  is  a  reason¬ 
able  compromise  and  it  leads  to  r„  —  877/V.  However,  as  we 
shall  see,  to  increase  the  resolution  we  might  use  a  larger 
value  for  f„. 

The  accuracy  of  the  method  and  the  attainable  resolution 
depend  on  the  relationship  between  the  unknown  frequencies 
cjf  and  the  sampling  frequency  u>0.  If  all  unknown  frequencies 
are  multiples  of 

w,  *  f/UJ0 

then  the  problem  is  essentially  discrete.  If  the  unknown  fre¬ 
quencies  and  their  differences  are  large  compared  to  co0,  then 
the  error  is  small  because  it  is  of  the  order  of  ui0. 

The  problem  of  determining  u>(  is  difficult  if  w0  is  of  the 
order  of  w(,  and  uif  is  not  an  integer  multiple  of  ui0 

W(a(r,  +  a)ui„  Icxl  <  -j . 

In  this  case,  the  resolution  error  u>„/2  is  of  the  order  of  u,t. 
Furthermore,  aliasing  generates  spurious  frequencies  in  the 
vicinity  of  uif.  Indeed,  if 

/(»)-«,“'f 


then 

W  =  gl2ir/N 

yielding  the  discrete  spectrum  (Fig.  5) 


.v-i 

Fm  =*  £  fn 

n*  0 


(m-r;-a)/V 
1  -  W  ‘ 

l-w(m-r'-a» 


To  improve  the  accuracy,  we  can  repeat  the  process  with  a 
larger  value  of  t„,  using  as  starting  B0  the  set  containing  only 
the  estimated  frequencies  u>,  and  their  neighbors. 

IV.  Illustrations 

We  illustrate  the  method  with  several  examples  involving 
signals  whose  unknown  frequencies  cannot  be  determined 
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Fig.  6.  Given  segment  wx  (r)  and  its  Fourier  transform  Wx(uj). 


Fig.  7.  Result  of  the  iteration  for  n  »  20  and  n  »  70. 
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Fig.  8.  Chen  data  segment  for  S/N  »  15, 11,  and  5  dB. 


from  direct  Fourier  analysis.  In  these  illustrations  we  con¬ 
sider  several  noise  levels.  With 

>vi(f)a/(f)  +  /»(r) 

the  given  data,  we  define  the  signal-to-noise  ratio  S/N  as  the 
ratio  of  the  energies  of  f(t)  and  n(r)  in  the  data  interval.  In 
ail  examples,  the  noise  is  white  and  is  uniformly  distributed 
:n  the  interval  (-c  to  c).  The  ratio  S/N  is  changed  by  changing 
the  size  of  c. 

The  computations  are  earned  out  with 

;V*256  /„  =  1  Hz  ta  *  1/256  s. 

To  avoid  large  scaling  factors,  we  divided  all  frequency  com¬ 
ponents  by  jV/2.  In  the  examples  we  show  also  the  value  of 
the  parameter  p  [see  (16)]  and  of  the  initial  threshold  level 

<i- 

Example  1:  The  unknown  signal  is  a  sum  of  two  sine  waves 


In  Fig.  6  we  show  the  given  segment  of  the  unknown  signal 
and  its  spectrum.  As  we  see  from  the  figure,  the  frequencies 
/,  and  fi  are  not  visible.  The  initial  threshold  is  e,  =>  0.15  and 
its  value  at  the  nth  iteration  is  obtained  from  (16)  with  p  = 
0.99.  In  Fig.  7  we  show  the  results  of  the  iteration  for  n  ■  20 
and  n  =  70.  At  the  70th  iteration  the  frequencies,  amplitudes, 
and  phases  of f(t)  are  recovered  exactly. 

We  note  that,  in  this  case,  the  values  of  e,  and  p  are  not 
critical.  Any  value  of  p  between  0.9  and  0.99  and  of  e,  be¬ 
tween  0.05  and  0.15  is  adequate.  The  iteration  was  per¬ 
formed  also  with  a  data  interval  containing  M  *  4 1  sampling 
points.  In  this  case,  the  results  are  sundar  but  the  speed  of 
convergence  is  slower. 

b)  We  consider,  next,  noisy  data  with  various  S/N  ratios  as 
in  Fig.  8.  in  all  cases, 

M. 3  5 1  q  3  0.99  e,  ”  0.1 5. 


fit)  *  1 .5  cos  (30irf  +  60°)  +  l  .25  cos  (20irr  +  30°) 

and  the  unknown  frequencies  /,  *  10  Hz  and  /j  =  IS  Hz  are 
integral  multiples  of  the  sampling  frequency 
a)  We  first  assume  that  the  data  interval  contains  M  -  5 1 
sampling  points  and  rr(r)  =  0. 


The  iteration  was  performed  several  times  with  the  same 
signal  but  with  different  samples  of  noise.  As  the  following 
indicates,  the  results  are  not  the  same  for  all  samples:  S/N  a 
15  dB  (c  =  0.375).  Six  samples  were  tried.  In  five  of  these, 
the  frequencies/,  and /2  were  found  exactly.  S/N  *  11  dB 
(c  =  0.625).  Fourteen  samples  were  tned.  In  nine,  we  ob- 
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U)  (b) 

Fig.  9.  (a)  Given  segment  w,(r)  of  /(f).  (b)  Fourier  transform  of  w,(rj. 


(a)  (b) 

Fig.  10.  Result  of  the  iteration  for  n  »  30  and  n  *  100. 


tuned  /,  and  /2  exactly.  In  four  cases,  an  error  of  1  Hz  de¬ 
veloped.  In  one  case,  the  iteration  yielded  not  two  but  three 
frequencies:  /,  =*  9  Hz,  /2  *  14  Hz,  and  /3  *  15  Hz.  S/yV  =  5 
dB  (c  *  1 .25).  This  is  a  very  noisy  case.  Of  the  eleven  samples 
tried,  three  gave  the  correct  answer,  two  yielded  1  Hz  error, 
five  resulted  in  2  Hz  error,  and  in  one  case  the  frequency  /2  ■ 
IS  Hz  was  lost. 

Example  2:  In  this  example  /(f)  consists  of  three  sine  waves 
and  the  data  are  noiseless.  We  consider  two  cases.  In  the  first 
case,  the  unknown  frequencies  are  multiples  of  oj0.  In  the 
second  case,  they  are  not. 

a) 

/(f)  ■  1.5  cos4irf+  1.5  cos(18irf +  60°) 

+  1.25  cos  (28irr  +  30°). 

We  start  with  the  following  values  of  the  relevant  parameters: 

M  3  59  p  =*  0.95  e,  =»  0.20. 

In  Fig.  9  we  plot  the  given  segment  /(f)  and  its  spectrum. 
Fig.  10  shows  the  results  of  the  iteration  for  n  =  30  and  n  = 
100.  At  the  1 00th  iteration  the  frequencies,  amplitudes, 
and  phases  of  /(f)  are  recovered  exactly . 

Again  the  values  of  p  and  «i  are  not  critical.  Essentially 
the  same  results  are  obtained  if  the  data  interval  is  reduced 
to  M  *  5 1  provided  that  p  is  not  less  than  0.95 . 

The  method  has  been  tried  also  for  a  smaller  data  interval. 
However,  the  convergence  is  slow  and  the  result  inaccurate. 
With  .Vf  =  4l,  p  *  0.99,  *0.20  the  component  with  the 
lowest  frequency  is  lost. 

b) 

/ff)«  1.5  cos  4.8irf  +  1.5  cos(18irr  +  60°) 

+  1.25  cos  (29. 2f  +  30s). 


In  this  case. 

/i  -  (2  +  0.4) /„  ft  *  9 f„  h  *  (14  *  0.6) /„. 

We  usedM  =  59, m  *  0.95,  and  e,  *  0.20. 

With  an  FFT  size  yV  *  256,  we  obtained  after  350  iteration 
steps  the  frequencies  2  Hz,  9  Hz.  and  15  Hz  (Fig.  1  lc). 

Increasing  the  FFT  size  to  .V  =  512,  we  found  in  200  steps 
the  frequencies  2.5  Hz,  8.75  Hz.  9.25  Hz,  and  14.5  Hz.  (Fig. 
lid). 

We  note  that  the  accuracy  in  the  evaluation  of  coefficients 
of  different  levels  can  be  improved  if  the  threshold  level  e„ 
is  not  constant  through  the  band  but  it  takes  different  values 
in  the  vicinity  of  each  frequency.  This  is  demonstrated  in 
the  next  example. 

Example  3:  The  unknown  signal  is  a  sum  of  five  sine  waves. 
/(f)  *  1 .5  cos  4irf  +  1 .25  cos  ( 1 2irf  +  30s) 

+  0.375  cos  (40irf  +  60s)  *  0.625  cos  50jrf 
+  1 .25  cos  (60irr  +  45s) 

with  frequencies  2,  6.  20.  25,  and  30  Hz;  the  noise  is  zero. 
In  Fig.  1 2  we  show  the  given  data,  obtained  with  M  -  71 .  and 
their  spectrum.  In  the  iteration  we  assume  that  p  *  0.99  and 
e,  =  0.04.  The  level  of  the  threshold  level  at  the  nth  iteration 
is  defined  as  in  (16).  However,  it  is  not  constant  throughout 
the  band.  Its  value  is  determined  from  the  behavior  W„(w) 
in  the  vicinity  of  each  maximum  (Fig.  13). 

In  Fig.  13  we  show  the  iteration  for  n  *  10  and  n  *  20.  At 
the  50th  step,  (Fig.  14)  we  recover  the  frequencies  2.  6.  25. 
and  30.  As  it  is  clear  from  the  figure.  Wn( a>)  contains  a  peak 
in  the  vicinity  of  /■  20.  To  determine  its  exact  location  we 
introduce  the  following  variation  to  the  method:  we  subtract 
from  the  given  data  the  recovered  portion  of  /(r)  and  repeat 
the  iteration  starting  with  the  new  data  J(t )  so  obtained. 
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Fig.  11.  (a)  Given  segment  W|(f).  lb)  Founer  transform  of  w,  (r). 
(c)  Result  of  the  iteration  for  n  *  350  and  FFT  size  At  »  256.  (d)  Re¬ 
sult  of  the  iteration  for  it  *  200  and  FFT  size  M  •  5 12. 

In  Fig.  IS  we  show  d(t)  and  its  spectrum  Z>(ui).  The  unknown 
frequency /=  20  is  recovered  at  the  20th  step  (Fig.  16). 

The  iteration  was  performed  also  with  a  smaller  data  seg¬ 
ment  (Af  =  61).  The  results,  however,  were  similar  but  the 
convergence  slower. 

Example  4:  To  test  the  limits  of  the  method,  we  consider 
as  a  last  case  an  example  where  the  data  interval  is  less  than 
one-half  the  unknown  period,  and  the  unknown  frequency  is 


IIMtl 


(b) 

Fig.  12.  (a)  Given  segment  u^lr).  (b)  Founer  transform  of  W]  <r). 


0  20  — —  «0 

1 1  Hal 


Fig.  13.  Result  of  the  iteration  for  n  *  10  and  it  ■  20. 

not  a  multiple  of  ui0  so  that  the  aliasing  is  significant.  We  as¬ 
sume  that 

/(f)  =  1 .25  cos  (5  Ant  +  30°)  T  =  0.08  s. 

This  yields  M  =  41  sampling  points  in  the  data  interval. 

The  iteration  was  performed  with  p  =  0.99  and  ei  =0.05. 
We  considered  four  different  signal  levels  (Fig.  IT), 
a)  n(t)  =  0.  At  the  40th  iteration  we  recover  the  frequency 
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0  *0  - -  40 

f  (Hi) 

Fig.  14.  Reailt  of  the  iteration  for  n  *  50. 
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Fig.  15.  (a)  New  data  segment  d(t ).  (b)  Fourier  transform  of  dit). 
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((Hi) 

Fig.  16.  Result  of  the  iteration  for  n  »  20. 


Fig.  17.  Gwen  data  segment  wt(f)  for  various  noise  levels,  (a)  SIN  - 
22  dB.  (b)S/N*  12  dB.  (c)  S/N  •  8  dB.  (d)  SIN  •  2  dB. 


/”  3  Hz.  Thu  is  the  nearest  sampling  frequency  to  the  un-  quency  and  its  two  neighbors  /=  2.5  and  f*  3.5.  After  n  = 
known  f\  *  2.7.  However,  since  the  resolution  frequency  150  steps,  we  recovered  the  frequency  /•  2.5  (nearest  to  the 
/,J1  is  of  the  order  of  flt  the  error  is  large.  To  reduce  it,  unknown  f\  *  2.7). 

we  increase  the  sampling  interval  from  ta  =  1/256  to  ta  -  The  process  was  repeated  with  ta  *  1/64,  that  is,  for  f0  = 
1/128.  This  yields  fa  =  1/2  Hz  but  the  number  of  sampling  1/4  and  .Vf  *  11  sampling  points.  The  iteration  yielded  two 
points  is  reduced  to  Af*21.  The  iteration  starts  from  the  frequencies:  /„  =  2.5  and  fb  =  2.75  with  amplitudes  |ca|  = 
band  B0  consisting  of  the  location  /*  3  of  the  recovered  fre-  0.616,  |c4|  »  0.630  and  phases  *  30.06°.  =  31.44°,  re- 
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TABLE  1 

Results  of  the  Iteeation  foe  20  Noise  Samples 
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Fig.  18.  wi(f)  »  1.25  cos  (5.4»f  +  30°)  ♦  n(t)  |r|  <  0.08  s.  (a)  Fourier  transform  of  W|<f).  <b>  Result  of  the  iteration 
for *  1  Hz  (i»f»  41)  and  n  *  100.  /(r) «  1.38  cos  (6irr  +  25-5°).  (c)  Result  of  the  iteration  for/„  >0.5  Hz  (M  *  21) 
and  n  »  200.  /(f)  *  1.23  cos  (5wf  +  3 1.8°).  (d)  Result  of  the  iteration  for  /„  »  0.25  Hz  (A#  «  11)  and  n  *  200.  /(f)  ■ 
0.62  cos(5»f+  30.1®) +  0.63  cos(5.5trf  +  31.4°)  =.  1.246  cos  (2.63trf  +  30.7°). 


spectively.  The  location  /  and  amplitude  c  of  the  unknown  In  Fig.  18  we  show  the  results  for  S/N  3  8  dB  and  f0  =  1. 
frequencies  was  finally  estimated  by  interpolation,  yielding  0.5, 0.25  Hz. 


,  l<*»l 

lcal  +  |c»| 


=  2.63  Hz 


c  *ca  *c„  =  1.246L30.76®. 


b)  fi(r)#0.  We  considered  four  different  signal-to-noise 
ratios.  All  cases  were  performed  20  times  using  different 
samples  for  the  noise.  The  statistical  conclusions  are  de¬ 
scribed  in  Table  I.  The  numbers  in  the  table  are  the  esti¬ 
mates  m  Hz  of  the  unknown  frequency.  Numbers  in  paren¬ 
theses  indicate  in  how  many  samples  out  of  20  the  cited 
estimates  were  obtained. 
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Theory  and  Applications  of  Running  Transforms 
6.  1.  Introduction 

In  the  early  years  of  filter  design,  the  objective  was  to  find  a 
realizable  system  whose  frequency  response  H(uu)  met  certain  requirements 
(low-pass,  high-pass,  or  band-pass,  for  example).  The  input  f(t)  was  viewed 
as  a  deterministic  signal  characterized  in  terms  of  its  Fourier  transform 
F(x)  and  the  output  g(t)  was  evaluated  from  the  inversion  formula  [LI 


g(t)  = 


1 


F(iu)H(ui  )ejujt  duu 


(1) 


This  is  still  a  design  method,  however,  in  the  last  two  or  three  decades, 
statistical  considerations  (optimum  estimation  or  detection,  for  example)  led 
often  to  design  requirements  based  on  the  formula  [2] 


EC  g2(t)}  =  J  S(uu)|H(tij)|2  dtu  (2) 

-00 

expressing  the  average  intensity  of  the  output  g(t)  of  the  filter  in  terms  of 
the  power  spectrum  S('ju)  of  the  input  f(t). 

In  both  cases,  the  design  involves  frequency  domain  considerations. 
However,  the  pertinent  spectra  are  global:  In  (1)  F(uu)  is  the  Fourier  transform 
of  the  entire  signal  f(t);  in  (2)  it  is  assumed  that  f(t)  is  stationary,  i.  e.  ,  its 
statistical  properties  are  the  same  for  all  t. 

In  a  number  of  applications,  it  is  desirable  to  design  filters  whose 
parameters  depend  not  on  global  but  on  local  properties  of  the  signals  to  be 
processed.  This  leads  to  a  time- varying  system 

adaptively  controlled  according  to  some  criterion  based  on  local  properties  of 
the  input. 


1 
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A  typical  non-recur sive  discrete  realization  of  such  a  system  is 
shown  in  Fig.  1.  The  input  x[n]  is  a  sequence  whose  N  values  x[n],  . .  .  , 
x[n-N+l]  form  the  input  vector  X[n]  .  The  resulting  output  g [n!  is  a 
weighted  sum 

N-l 

y[n]  =  )  a^[n]  x  [n- k]  =  XT[n]  A[n]  (3) 

k=0 

where  A[n]  is  a  vector  whose  N  components  are  the  time-varying  weights 
a^Tn]  of  the  system. 

The  above  system  has  finite  memory,  hence,  its  response  y[n] 
depends  only  on  local  properties  of  the  input  x[n].  However,  since  it  is 
essentially  a  time-domain  filter,  it  does  not  utilize  directly  the  frequency- 
domain  properties  of  x[n]  .  For  a  number  of  applications  (noise  or  clutter 
rejection,  for  example)  this  is  a  disadvantage.  Another  disadvantage  is  the 
fact  that  the  number  N  of  its  adaptively  controlled  parameters  equals  the 
length  of  its  memory.  This  number  can  be  reduced  if  the  processing  is 
performed  recursively.  However,  adaptively  controlled  recursive  filters 
have  other  disadvantages. 

In  this  paper,  we  show  that  localized  frequency  domain  processing  is 
possible  if  the  spectral  properties  of  the  input  f(t)  are  expressed  in  terms  of 
its  running  transform  F(t,  jj).  The  main  objective  of  the  paper  is  the 
development  of  running  transforms.  However,  to  demonstrate  their  value  in 
filter  design,  we  include  two  applications. 

The  results  are  presented  in  analog  form.  Discrete  and  two-dimensional 
signals  are  considered  only  briefly. 


2 


2.  The  running  Fourier  transform 


c 

Fc(t,uu)  =  J  f(t+T)  e"^T  dT  (4^ 

-c 

where  c  is  a  given  constant.  For  a  fixed  t,  Fc(t,  x)  is  the  Four.er 
transform  in  the  variable  r  of  the  segment  f(t+T)  of  f(t)  shown  in  Fig.  2. 

The  concept  of  the  running  transform  is  not  new  [4],  However,  in 
earlier  years,  it  was  of  limited  use  because  of  the  then  underlying 
computational  difficulties. 

The  function  Fc(t,uu)  describes  the  local  spectral  properties  of  f(t)  and 
their  variation  with  time.  As  c  increases  from  0  to  m ,  it  changes  from  a 
time  signal  to  a  frequency  signal.  At  the  two  extremes, 

F,  (t,  jo>  =  F(w)  ±  F Jt,  uu)  _  f(t)  (5) 

c  — 0 

We  shall  use  the  notation 

f(t)< - ♦  Fc(t,  jj) 

to  indicate  that  Fc(t,  jj)  is  the  running  transform  of  f(t).  The  subscript  c 
will  atten  be  omitted. 

From  the  definition  it  follows  readily  that 

f(t-tQ)  i - 1  F(t-to,  ju)  (6) 

jui.t  ju)  t 

e  c  f(t)  e  F(t,UHJJc) 
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A  3imple  illustration  follows: 


2  sine 


UU 


uu 


jouct  juu  t  2sinc(uu-uu  ) 

e  | - 1  e - S— 


UU-'JU 


More  generally,  if  f(t)  is  a  periodic  function  with  period  2c: 


Ar  = 


jr!uot 


_  JT_ 
1  o  c 


r=-. 


then 


2  sin  c(x  -r  ju  )  jruu  t 
F(t,uu)  =  ^  Ar  - — - 2.  e  0 


jd—  r  jo 


(7) 


r=-« 


Hence  [see  (6)]  , 


i*c‘  X  .  irV 

*  l  Ar  = 

r=-«» 


juu  t  -  2  sin(uu-uu  -rx  )  jruu  t 

e  'j  A  - E - 0  (8) 


T--m 


r  ju-uu  -ruu 
c  o 


This  result  will  be  used  later.  We  note  that  for  the  example  in  (8), 


j(uj  +  nuj  )t 

F(t,  uu  +  nuj  )  =  2c  A  e 
c  o  n 


(9) 


Properties.  Suppose  first  that  f (t )  is  a  deterministic  signal  with  Fourier 
transform  F(w).  In  this  case,  F(t,  uu)  has  a  Fourier  transform  in  the  variable 
t: 

m 

T“(w,  ju)  =  J  F(t,  jj)  e"^wt  dt  (10) 
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From  (4)  and  (10)  it  follows  readily  that 


c 

Ttw.u,)  =  f  F(w)ejWr  e'jWTdr=  Ffw) 

J  W-(JU 

-c 

This  shows  that  for  a  specified  yu ,  F(t,  yu)  can  be  considered  as  the 
output  of  a  linear  system  with  input  f(t)  and  system  function 

j(w>ttl)  = 

W-U) 

where  w  is  the  frequency  variable  and  u>  a  constant  parameter  (Fig.  3a) 

We  now  assume  that  f(t)  is  a  stationary  random  process  with  auto¬ 
correlation  R(t)  and  power  spectrum  S(w).  In  this  case,  F(t,ui)  is  also  a 
random  process  stationary  in  the  variable  t  and  non- stationary  in  the 
variable  ou .  We  shall  determine  its  autocorrelation 

Rpp(r,ij:)  =  E{F(t+T,uu)  F  (t.uu)} 
and  power  spectrum 

OB 

SFF^*VJ^  =  I  tu)e'^"  dr 

-  m 

Since  F(t,  jj)  is  the  output  of  the  system  of  Fig.  3a  with  input  f(t), 
it  follows  from  a  well  known  theorem  H2]  that 


S-,r(w,  jj)  =  S(w ) 


4  sin  c(w-x) 


The  inverse  transform  of  the  above  in  the  variable  w  yields 


2c 

=  j  (2c- jx|)  R(x+T)e"^wxdx 
-2c 

This  follows  from  the  convolution  theorem  and  the  fact  that  the 
inverse  transform  of  the  fraction  in  (15)  is  a  triangular  pulse  [1] 
multiplied  by  e"-1"*. 

For  a  fixed  t,  F(t,uj)  is  a  non- stationary  process  in  uu .  We  shall 
show  that 


E(F,t,  -  f  J  S(w)  1  dw 


Clearly,  F(t,  u)  and  F(t,  v)  are  the  outputs  of  two  systems  with 
common  input  f(t)  and  system  functions 


2  sinc(w-u) 
w-  u 


and 


2  sin  c(w-v) 
w- v 


respectively  (Fig.  3b).  Hence,  their  cross-power  spectrum  in  the  variable 
t  equals 

,  4  sinc(w-u)sinc(w-v) 
b(W) - [w"-u)(w-v) - 


This  leads  to  the  conclusion  that 


E[F(t+r.  u)F::!(t,  V)}  =  J  S(w)  ejwT  dw 


(w-  u)(w- v) 


and  (17)  follows  with  t=0  . 
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We  note  finally  with  u=v=uu  that 


E  (  I  F(t,  „)|2]  =  i-  f  S(w)  ■=(»-»)  d„  (1„ 

(w-uu) 

formula.  We  shall  now  show  that  f(t)  can  be  expressed  in  terms 
of  the  samples  F(t,  mcu^)  of  its  running  transform: 

00 

f(t)  =~2F  1  F(t,  m(uo)  (20) 

m=-» 

Proof.  We  expand  the  function  f(t +r)  into  a  Fourier  series  in  the  variable 
t  •  The  coefficients  of  this  expansion  in  the  interval  (-c,c)  are  given  by 
[see  (4)] 


i  i*C  r  , 

TF  J  f(t+T)e  dr  =  JE  F(t*m‘U0) 

-c 


(21) 


1  jmtu  r 

f(t+i-)  = 2  F(t,majo)e 

m=-» 


(22) 


for  |  t |  <  c  .  And  with  r=0,  (20)  follows. 

Sampling  expansion.  The  preceding  result  is  a  version  of  the  sampling 
theorem  and  is  based  on  the  fact  that  F(t,jj)  is  the  Fourier  transform  of 
the  time-limited  function  shown  in  Fig.  2.  Thus,  f(t)  and,  therefore, 
F(t,  (i) )  is  uniquely  determined  in  terms  of  F(t,  maj0).  In  fact,  F(t,  x)  can 
be  expressed  explicitly  in  terms  of  its  samples  F(t,  m«o): 


F(t,  uu)  =  ) 


F(t,  m(jjo) 


m=  -« 


sin  c(a)-mcto) 
c('ju-mcuQ) 


(23) 
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_Proof.  Inserting  (22)  into  (4)  and  integrating  term-wise,  we  obtain  (23). 

Interpolation.  From  (23)  it  follows  with  uu=  r<u  +  m  / 2  that 
-  o  o 

"  r_k 

F(t,r.x0+  xq/2)=  i-  V  ^'Vl/2  <24) 

k— — to 

This  result  can  be  used  to  increase  the  resolution  in  the  spectral 
representation  of  f(t).  Indeed,  suppose  that  we  know  the  samples  Fc(t,myjQ) 
of  Fc(t,  iju).  If  we  double  the  interval  (-c,c),  we  obtain  the  running  transform 
F2c(t,  uu )  whose  samples  are  F2c(t,  mu)  Q/2).  We  shall  express  these  samples 
in  terms  of  Fc(t,mujo). 

From  the  definition  (4)  it  follows  readily  that 

F2c(t.w)  =  Fc  (t-c,u>)  +  Fc(t+c,0B)  (25) 

If  m  =2r  is  even,  then  (2  5)  yields 

F2c(t,muuQ/2)  =  Fc(t-c,ru,o)  +  F^t+c.ru^)  (26) 

If  m  =  2r  +1  is  odd,  then  [see  (24)] 

®  r-k 

F2c(t,  m„0/2)  •  —  [Fc(t-c,  k%)  V  Fc(t+c,  ta»0)]  (27) 

k=  — 30 

Reasoning  similarly,  we  can  express  the  samples  F  (t,  mouo/n)  of 
Fnc(t,uu)  in  terms  of  the  samples  Fc(t,m(juQ)  of  Fc(t,au)  . 
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Existance  and  uniqueness.  An  arbitrary  function  F(t,  jj)  of  two  variables 
is  not,  in  general,  a  running  transform.  It  must  satisfy  the  sampling 
expansion  (23).  We  shall  examine  the  properties  of  its  samples 

am(t)  =  F(t*m“J0)  (28) 

Suppose,  first,  that  f(t)  is  a  deterministic  signal  with  Fourier 
transform  F(w).  From  (11)  it  follows  that,  in  this  case,  the  Fourier 
transform  Am(w)  of  am(t)  is  given  by 

2  sin  c(w-moi  ) 

Am(w)  =  F(w,m«o)  =  F(w)  -  - - -  (29) 

This  leads  to  the  conclusion  that 

A  (w)  =0  for  w=r r  4  m  (30) 

rn  o 

If  f(t)  is  a  stationary  stochastic  process  with  power  spectrum  S(w), 
then  arn(t)  is  also  a  stationary  process  and  its  power  spectrum  S^fw)  is 
given  by  [see  (15)] 

4  sin^  c(w-muu  ) 

sm(w)  =  SFF(w,m(jjo)  =  S(w)  - 2~°-  (31) 

(w-mujo) 

From  this  it  follows  that 

S  (w)  =  0  for  w=ruu  ri  m  (32) 

m  o  '  1 
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We  consider  now  the  converse  problem:  Given  a  function  a  (t) 

m  ' 

we  wish  to  find  a  function  f(t)  such  that  the  mth  sample  F(tfmtuQ)  of 
its  running  transform  equals  am(t).  We  assume,  of  course,  that  if 
a  (t)  is  a  deterministic  signal,  then  it  Fourier  transform  A  (w) 
satisfies  (30);  if  it  is  a  random  process;  then  its  power  spectrum  S  (w) 
satisfies  (32). 

To  determine  f(t),  we  form  the  inverse 


^(w,  mu)Q)  = 


.  sm  c  (w- moj 


of  the  system  J(w,muuQ)  defined  in  (12).  We  maintain  that  if  the  input  to 
this  system  is  the  given  function  am(t),  then  the  resulting  output  is  the 
unknown  f(t). 

Indeed,  because  of  our  assumptions,  f(t)  exists  [see  also  (11)  and  (15)]  . 

Furthermore,  the  mth  sample  F(t,mujo)  of  its  running  transform  is  the 

output  of  the  filter  J (w,  m<uQ)  with  input  f(t),  that  is,  it  is  the  output  of  the 

filters  ^(w.muu  )  and  J(w,  muj  )  is  cascade  with  input  a  (t).  Hence, 
oo  r  m 


F(t,  mm  )  =  a  (t) 
o  m 


We  conclude  with  the  investigation  of  the  uniqueness  of  the  signal 
f(t)  30  obtained.  From  this  discussion  it  will  follow  that  the  functions 
F(t,mujo)  are  not  arbitrary. 

Suppose,  first,  that  the  running  transform  F(t,  *)  vanishes  identity  in 

t  for  some  ju  =  .u  : 

c 


F(t,  juc)  ■  0 


(34) 
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If  f(t)  is  a  deterministic  signal  with  Fourier  transform  F(w),  then 
(34)  yields  [see  (11)] 


2sinc(w-'ju  ) 

)  =  F(w)  — — - — 

c  w-uu„ 


(35) 


This  is  an  identity  in  w  and  can  be  satisfied  only  if 


F(w)  bf  4(w-(uc-ruuo) 


(36) 


where  the  prime  in  the  summation  indicates  that  r  takes  all  integer 
values  except  r=0.  From  (36)  it  follows  that 


f(t)  = 


1  JV  v  '  . 

■Zff  6  1  br  6 

r 


(37) 


Consider,  next,  two  signals  f^(t)  and  f2(t)  whose  running  transforms 


are  identical  in  t  for  some  oj  =  jj  : 

c 


Fx(t.  tuc)  *  F2(t,  uuc ) 

With 

f(t)  =  fx(t)-f2(t)  l - 1  F  1(t,  uu)  -  F2(ttiii)  =  F(t,w) 

it  follows  that  F(ttU)c)  ■  0.  Hence  [see  (37)]  , 

l  jot)  t  '  jrou  t 
^(t)  =  f2(t)  +  ^  e  )  br  e  ° 


(38) 


(39) 
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This  shows  that  if  the  running  transform  F(t,  tu)  of  f(t)  is  specified 
for  some  x  =  myj  as  in  (28),  then  f (t )  is  unique  within  a  periodic  function 
of  period  2c  whose  mth  harmonic  is  zero.  Thus,  the  running  transforms 
of  the  functions 


i  jrsju  t 

«‘>  +  Tn  V  br*  » 

r^rn 


given  by  [see  (8)] 


l  br 
r^m 


sinc(uu.-riju0)  jrx  t 
,J)  -rsuo  6 


(41) 


are  identical  in  t  for  uu  =moj  . 

o 

Similar  conclusions  hold  if  f(t)  is  a  random  process. 

3.  Filtering 

The  applications  of  running  transforms  in  filtering  and  estimation 
are  based  on  the  following: 

If  a  linear  system  is  time-limited,  i.  e. ,  if  its  impulse  response 
h(t)  is  such  that 


h(t)  =0  for  ]t|  >  c 


(42) 


then  the  response 

c 

g(t)  =  |  f(t-o)h(a)da  (43) 

-c 


12 


-60 


to  an  arbitrary  input  f(t)  is  given  by 


*<*>  -  h  1 


F(t,  mouo)  H(muuo) 


m=- 1 


where 


(44) 


c 

H(tu)  =  J  h(t)e'ju,tdt  (45) 

-c 

is  the  system  function  and  F(t,uu)  is  the  running  transform  of  the  input. 

Proof.  Setting  r  =  -a  in  (22),  inserting  into  (43)  and  integrating  term-wise, 
we  obtain  (44). 

Equation  (44)  expresses  g(t)  in  terms  of  the  samples  H(muuo)  of  H(uj) 
and  the  spectral  properties  of  f(t)  in  the  interval  (-c,  c).  It  involves  an 
infinite  number  of  terms,  however,  only  a  small  number  of  these  terms 
need  be  considered.  Indeed,  if  f(t)  is  a  deterministic  signal,  then 


F(t,  muuQ)  = 


F(w) 


sin  c(w-muuo) 
rr(w-muu0) 


dw 


(46) 


This  follows  from  (11)  and  the  Fourier  inversion  formula.  If  f(t)  is  a 
random  signal,  then  [see  (19)1 


-  sin  c(w-muo  ) 

E(  |F(t,m»a)n«  \  S(w) - dw  (47) 

(w-muj0) 


Both  relationships  show  that,  if  the  bandwidth  of  f(t)  equals  3,  then 

the  number  of  significant  terms  F(t,  muuo)  in  (44)  is  of  the  order  of  j/’JUq  . 

This  number  can  be  made  as  small  as  desired  if  the  constant  c  =rr/'JU  is 

o 

sufficiently  small. 
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Running  low-pass  and  band-pass  filters.  A  special  case  of  (44)  is  the 
low-pass  filter  obtained  with 


H(mui0)  - 


1  |m  [  <  M 

0  |m  I  >  M 


(48) 


In  this  case,  (44)  yields 
M 


M  M 

=  7c  Y  mujo^  =  Tc  Ft't‘  +  4"^  ReF(t,muj0) 


(49) 


m=-M 

Inserting  (4)  into  (49) ,  we  obtain 
M  c 


m=l 


.  .  1  {  r°  jmit  a  c  M  jmau  a 

g(t)  =  2c  \  j  f(t"a)e  da=  J  f(t-a)  ^  e  0  da 


m=-M  -c 


■c  m  =  -M 


This  shows  that  the  mpulse  response  of  the  running  low-pass  filter 
is  given  by 


M 

,  *_  jmui  t  sin(M+l/2)uu  t 

•  --z-.nl.  t/l 

m=-M  ° 


for  |t  j  <  c  and  it  equals  zero  for  It  |  >c  . 

The  corresponding  frequency  response  is  given  by 

M 


(50) 


sinc(uu-mu)  ) 

H-.fa)  =  )  - , - ° 

™  -  cfou-mcu  ) 

m=-M  ° 


(51) 
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In  Fig.  4,  we  plot  the  functions  ’a^j(t)  and  H,^(u))  for  a  fixed  c  and 

for  M  from  0  to  five.  The  impulse  response  h^(t)  is  zero  for  |t)  >  c; 

however,  its  effective  duration  is  smaller.  Using  as  its  measure  the 

width  t  of  the  main  lobe,  we  find 
c 


,  _  2tt  _  2c 

c  (2M+1)(juq  2Mt1 


(52) 


The  frequency  response  ti^('oj)  is  a  low-pass  curve  with  cut-off 
frequency 


'JU 


c 


(53) 


In  Fig.  5,  we  plot  H^(t )  and  H^(uu)  for  a  fixed  value  of  tc  and  for 

M  from  0  to  five.  The  bandwidth  uu  of  the  filter  or  the  duration  t  of 

c  c 

its  impulse  response  can,  thus,  be  controlled  by  varying  c  or  M. 

For  M=0,  (49)  yields 

c 

g(t)  =  F(t,  0)  =  J*  f(t+T)dr  (54) 

-c 

This  is  the  moving  average  of  f(t)  and  it  has  the  following  familiar 
applications:  Suppose  that 


f(t)  =  y(t)  +  n(t)  (55) 

where  y(t)  is  a  signal  to  be  estimated  and  n(t)  is  noise.  The  component  of 
g(t)  due  to  the  noise  decreases  if  c  is  large,  however,  the  distortion  of 
y(t)  (bias)  due  to  smoothing  increases.  To  minimize  the  resulting 
mean- square  error,  we  must  choose  c  adaptively  [5]  taking  into 
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3.  Discrete  processing: 

The  numerical  implimentation  of  the  running  Fourier  filter  involves 
truncation  of  the  sum  in  (44)  and  evaluation  of  the  samples 

F[n,  m]  =  F(nT,mu)  ) 
o 

of  the  running  transform  F(t,  jj)  in  terms  of  the  samples 


f[n]  =  f(nT) 


(60) 


of  f(t).  In  the  above. 


T=  2c /N 


(61) 


is  the  sampling  interval  and  it  is  chosen  so  as  to  yield  a  small  aliasing 
error.  The  truncation  and  aliasing  errors  can  be  determined  a3  usual 
using  the  inversion  formula  [see  (11)  ] 


F(t,  jo) 


|  F(») 


sinc(w-'ju)  jwt  , 
-ff(w-ou)  e  dw 


(62) 


for  deterministic  signals  and  the  relationship  [see  (15)] 


EC  |F(tf!U)|2}=  f  S(w)dw 

rr(w-ctr 


(63) 


for  random  signals. 
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We  present,  next,  the  digital  version  of  running  transforms. 

The  running  discrete  Fourier  series  (abbreviation:  DFS)  of  a 
discrete  signal  f[n]  is  the  sum  [3]  ,  [6] 

N-l 

F[n,m]  =  )  f[n-k]wkm  w  =ej2ft/N  (64) 

0 

where  N  is  a  given  integer.  For  a  fixed  n,  F[n,  ml  is  the  DFS  in  the 

of  f[n] 

variable  k  of  the  segment  f[n-k]/shown  in  Fig.  7.  Unline  (4),  we  chose 
the  discrete  time  n  not  in  the  center  but  at  the  end  of  the  segment  shown 
in  Fig.  7.  We  did  so  to  conform  with  the  usual  notations  for  DFS.  If  (64) 
is  used  as  the  discrete  version  cf  (4),  a  delay  equal  to  c  =  NT/2  must  be 
allowed. 

Inversion  formula  .  From  the  inversion  formula  for  discrete  Fourier  series 
it  follows  that  [3] 

N-l 

f[n-k]  =  4-)  F[n,m]w’km  0<k<N  (65) 

m=0 

With  k=0,  this  yields 
N-l 

f[n]  =  }  F[n,m]  (66> 

m=  0 

Filtering.  Suppose  that  the  discrete  signal  f[n]  is  the  input  to  a  non-recursive 
filter  of  length  N.  We  maintain  that  the  resulting  output 
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N-l 

g[n]  =  S  f[n-k]  h[k]  (67) 

te=0 

is  given  by 

N-l 

g[n]  =  -rr  )  F[n,  m]  H(wm)  (68) 

m=0 

where 

N-l 

H(z)  =  ^  h[k]  z"k  (69) 

k=  0 

is  the  system  function  and  F[n,  ml  is  the  running  DFS  of  the  input. 

Proof.  Inserting  (65)  into  (67)  and  changing  the  order  of  summation, 
we  obtain 

N-l  N-l 

g[n]  =  )  P[n,  m]  )  h[k]  w_km  (70) 

m=  0  k=  0 

and  (68)  results  because  the  last  sum  equals  H(wm), 

The  filter  specified  by  (68)  has  been  used  in  a  different  context  as  a 
frequency  domain  interpolator  [7]  . 

Equation  (68)  is  the  discrete  version  of  (44)  and  it  forms  the  basis 
of  running  frequency  domain  filtering.  Its  implimentation  is  shown  in 
Fig.  8.  It  consists  of  a  memory-less  system  W  transforming  the  N 
samples  of  the  input  f[n]  into  their  DFS,  followed  by  N  weights  H(wm)  and 
an  adder.  The  discrete  version  of  the  running  low-pass  and  band-pass 
filters  discussed  earlier  results  if  Hlw31)  takes  the  values  zero  or  one. 
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The  system  of  Fig.  8  can  be  used  in  a  variety  of  applications. 

As  an  illustration,  we  present  briefly  the  frequency  domain  version 
of  the  Wiener  filter  stressing  the  advantages  over  the  familiar  approaches 
The  problem  under  consideration  is  the  estimation  of  a  discrete 
process  y[n]  in  terms  of  the  N  most  recent  samples  x[n-k]  of  another 
process  x[n]  (data).  The  usual  approach  to  this  problem  involves  the 
determination  of  the  N  weights  a^  such  that  if  the  sum 

N-l 

y[n]  =  )  a^  x[n-k]  =  [n]  A 

fc=0 

is  used  as  the  estimate  of  y[n]  ,  the  resulting  MS  error 
E  {  |yfrn  -  y[n]|2} 

is  minimum,  As,we  know  (orthogonality  principle),  the  optimum  vector 
A  is  given  by 


where 


R  =  E  {  X[n]  XT[nl  } 

is  the  correlation  matrix  of  the  data  vector  X[n]and 
T  -  Ef  X[n]  yCn]  } 
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is  the  cross-correlation  vector  between  the  data  vector  and  the  signal 
y[n]  to  be  estimated. 

This  approach  has  the  disadvantage  that  it  involves  the  inversion 
of  the  matrix  R  whose  order  N  is, in  general,  large.  This  problem  is 
particularly  severe,  if  R  is  not  known  in  advance  but  it  is  estimated  for 
each  n.  In  thi3  case,  R  depends  on  n(adaptive  Wiener  filter)  and  the 
inversion  must  be  performed  for  each  n. 

In  the  frequency  domain  version  of  the  Wiener  filter,  the  unknown 
signal  y[n]  i3  estimated  by  the  sum 

M 

TW-I  bm  F[n.m]  (73) 

m=-M 

where  F[n,  m]  is  the  running  DFS  of  the  data  sequence  x[n]  .  If  the 

optimality  criterion  is  the  minimization  of  the  MS  error,  then  the  vector 

B  of  the  unknown  weights  b  is  given  by 

m  a 

B  =  C"1  D  (74) 

where  C  is  a  matrix  of  order  2M+1  with  elements 

E  [  F[n,m.]  F*  [n,  m.]  } 
and  D  is  a  vector  with  elements 

E  f  F[n,  m]  y  [nj  ] 
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In  this  approach,  the  required  value  of  M  for  nearly  optimal 
estimation  is  small  compared  to  N.  Furthermore,  the  matrix  C  is 
nearly  diagonal  and  can  be  easily  inverted.  This  follows  from  the 
discrete  version  of  (17)  but  the  details  are  omitted. 

We  mention  without  elaboration  that  the  frequency  domain  Wiener 
filter  (73)  properly  modified  can  be  used  in  the  estimation  of  non- stationary 
processes. 

The  DFS  analyzer.  The  implimentation  of  the  frequency  domain  filter  of 
Fig.  8  involves  the  evaluation  of  the  running  DFS  F[n,m]  of  the  input.  It 
appears  from  (64)  that  the  required  number  of  multiplications  equal  N-l 
for  each  m.  However,  as  we  3how  next,  this  number  can  be  reduced  to  a 
single  multiplication  if  F[n,  m]  is  evaluated  recursively.  Indeed,  from  (64) 
it  follows  that 

F[n,  m]  =  wm  F[n-1,  m]  4  f[n]  -f[n-N]  (75) 

This  can  be  realized  with  a  first  order  system  if  we  have  access  to 
f[n]  and  f[n-N]  .  The  realization  of  the  component  W  of  the  running  filter 
of  Fig.  3  requires  2M+1  such  systems. 

4.  Two-dimensional  signals 

The  preceding  analysis  can  be  extended  to  two-dimensional  signals. 

We  conclude  with  a  brief  outline  of  typical  results,  omitting  all  proofs. 

The  running  Fourier  transform  of  a  signal  f(x,  y)  is  the  integral 
a  b 

F(x,  y;  u,v)  =  f  f  f(x+?,  y+ti)e"^u- +  vr|^  de  dri  (76) 

-a  -b 

where  a  and  b  are  given  constants. 
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Inversion  formula.  It  follows  as  in  (20),  that 

m 

f(x’  y)  =  TSZ  i  F<x’  V'  ml  V  rn2vo>  uo  =  a"  v0  =  T  <77> 
m1m2=  — 

Filtering.  If  a  two-dimensional  linear  system  is  space-limited,  i.  e.  , 
if  its  point  spread  h(x,  y)  is  such  that 


h(x,  y)  =0  for  |x|  >a  ,  |y|  >b 


(78) 


then  its  response 


a  b 


g(x,  y)  =  J  J  f(x-§,  y-r|)  d§  dn 


-a  -b 


to  an  arbitrary  input  f(x,  y)  is  given  by 


(79) 


S(x'  y)  =  TIB  1  F{x’  ^  mlUo'  m2vo}  H(ml  V  m2vo) 

in ^  ~ ~ ® 

where  H(u,  v)  is  the  system  function  and  F(x,  y;  u,  v)  is  the  running 
transform  of  the  input. 

Discrete  signals.  The  running  DFS  of  a  two-dimensional  discrete 
signal  f[n^,  r^]  is  the  sum 


Nrl  N2-l 


kjm  k2m2 


r  (jn^,  n2>  m^,  m21  =  ^  ^  ffn^-k^,  n2-k2'  W1  ”2 

k^=  0  k2=0 


(80) 


(81) 
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where 


jSfl/Nj  j2rt/N2 

w^  =  e  ,  w^  =  e 


Inversion  formula  .  It  follows  as  in  (65)  that 


Nrl  N2-l 


V  “23  8  to:  1  i  F[ni,n2;  m1,m2] 
^  “  _  /-\  _ 


m^=0  m2=0 


Filtering.  If  the  signal  ffn^,  r^]  is  the  input  to  a  non-recursive 
two-dimensional  filter,  then  the  resulting  output 


Nl-1  ^2_1 


g  L11!*  n2]  =  ^  ^  ftnj-kj,  n2-k2]  hO^,  k2] 

kj=0  k2=o 


is  given  by 


Nrl  N2-l 


1  c  r  m  m. 

g  ^2»  m2-  9  w2 


where 


V°  k2=o 


Nrl  N2-l 


—  -  *  k«  —  h 

H(zlf  z2)=  ^  )  hCkj,  k2]  z^  z2  2 

v°  V° 

is  the  system  function  and  F[n^f  n2;  m^,  m2]  is  the  running  DFS  of  the  input. 
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Section  7 

Adaptive  Frequency  Domain  Estimators 
7.1.  Adaptive  and  recursive  MS  estimators 

Adaptive  estimators  are  of  interest  because  they  require  no 
knowledge  of  prior  statistics  and  can  be  used  to  process  stationary  and 
non- stationary  signals  [1]  ,  [21  .  Recursion  is  desirable  because  it 
simplifies  the  required  arithmetic  operations.  In  this  paper,  we  develop 
a  new  estimator  that  combines  these  advantages  with  the  advantages  of 
frequency  domain  processing. 

To  place  our  development  into  familiar  perspective,  we  reexamine 
two  basic  methods  of  adaptive  filtering  using  as  illustration  the  estimation 
of  a  discrete  process  y[n]  in  terms  of  the  output 

N-l 

y[n]  =  ^  ak  [n]  x[n-kl  (1) 

k=  0 

of  a  non-recursive  time- varying  system  whose  input  is  the  data  process 
x[n]  .  In  vector  notation, 

y[n]=  XT[n]A[n]  (2) 

where  A[n]  is  the  vector  whose  components  are  the  N  weights  a^[n]  and 
T 

X  [n]  is  the  transpose  of  the  data  vector  X[n]  consisting  of  the  N  most 
recent  samples  x[n-k]  of  x[n]  .  The  resulting  estimation  error  is  given 
by 


e[n]  =y[n]-y[n]=  y[n]  -XT[n]  A[n] 


(3) 
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and  oar  objective  is  to  determine  the  vector  A[n]  so  as  to  minimize 
in  some  sense  this  error. 

The  adaptive  Wiener  filter.  If  the  processes  x[n]  and  y[n]  are  jointly 
stationary  and  the  optimality  criterion  is  the  minimization  of  the  MS 
value  of  e[n]  ,  then  the  optimum  A[n]  is  the  solution  of  the  system 
(orthogonality  principle  [3]) 

RA  =  r  (4) 

where  R  is  the  correlation  matrix  of  the  data  vector  and  p  is  the 
cross-correlation  vector  between  the  data  and  the  signal  y[n]  to  be 
estimated: 

R  =  E  [  X[n]  XT[n]  }  r=  E  {  X[n]  y[n]  ]  (5) 

The  solution  of  (4),  yields  the  familiar  Wiener  filter 

a  =  r'1  r  (6) 

This  filter  is  optimum  but  it  requires  prior  knowledge  of  the  matrix 
R  and  the  vector  p. 

For  the  estimation  of  a  process  with  unknown  statistics  we  can  use 
a  vector  A[n]  such  that 

R[n]  A[n]  =  p[n]  (7) 


3 


where 


n-1 

R[n]  =  (1-a)  )  z  X [n- k]  XT[n-k] 
k=l 


n-1 

f  [n]  =  (1-a)  ^  ak  X[n-k]  y  [n-k] 
k=l 

The  solution  of  (7) 

A[n]  =  R"1  [n]  r[n] 


0  <a <1 


(8) 


(9) 


is  the  adaptive  Wiener  filter. 

It  can  be  shown  (Appendix  A), that  under  general  conditions,  the 

A  a 

matrix  R[n]  and  the  vector  p  [n]  tend  to  R  and  p  respectively  as  n  -« 
and  a-  1.  Hence,  for  sufficiently  large  n  and  for  a  close  to  one,  the 
filter  specified  by  A[n]  is  nearly  optimum. 

We  note  that  the  adaptive  Wiener  filter  can  also  be  used  in  the 
estimation  of  non- stationary  processes  provided  that  their  statistics  vary 
sufficiently  slowly  (see  Appendix  A). 

For  a  given  n,  the  output  y[n]  of  the  above  filter  is  not  optimum  in 
the  mean- square  sense.  It  can  be  shown,  however,  that,  if  the  optimality 
criterior  is  the  minimization  not  of  the  mean- scuare  value  of  the  error  but 
the  minimization  of  its  weighted  time-average 
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n-1  n-1  ^ 

)  a*  e  [n-k]  =  )>  aK  I  y[n-k]  -  XT[n-k]  A[n-k]  ')  (10) 

r =i  k=i 

then  y[n]  is  optimum  for  any  n. 

The  realization  of  the  vector  A[n]  involves  the  evaluation  of  the 
vector  7  Ln]  and  the  matrix  R[n]  and  its  inverse  for  each  n.  The  underlying 
computations  can,  however,  be  simplified  if  they  are  carried  out 
recursively.  Indeed,  from  (8)  it  follows  readily  that 


R[n+1]  =  a R[n]  +°Va)  X[n]  XT[n] 

A 

f  [n+l]  =  ar  [n]  +^l-a)  X[n]  y[ni 

A 


(11) 


To  obtain  a  recursion  for  the  solution  A[n]  of  (7),  we  note  with  (11) 
and  (7)  that 


A  [n+l]  =  R-1  [n+l]  T  [n+ll  =  R_1[n+1]  {  af  [n]  Al-o.)  X[n]  y[n]  } 

A 

=  iR  '  [n+l]  R[n]  A[n]  -^(l-a  )  R[n+1]  X[n]  y  [nl 
=  R'1  [n+l]  |R[n+l]  -*(l-a)  X[n]XT[n]}  A[n] 

Al-a)  R_1[n+1]  X[n]  y[n] 
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y[n]  =  X^jh]  A[n] 


e[n]  =  y[n]  -y[n] 


we  conclude  that 


A[n+1]  =  A[n]  +  (l-a)  R_1[n+l]  e[nl  X[n]  (12) 

A 

We  can  similarly  derive  a  recursion  equation  for  the  inverse  matrix 
R_1[n+ll  . 

The  Widrow  filter.  A  simplified  form  of  (12)  is  the  algorithm 

A[n+1]=  A[n]  +  p.  e[n]  X  [n]  (13) 

where p  is  some  constant  that  determines  the  time  of  adaptation.  This  is 
the  basis  of  the  Widrow  filter  [4]  ,  [5]  and  it  is  justified  as  o.n  instantaneous 
version  of  gradient  seeking  least  MS  techniques  [5]  . 

Unlike  the  adaptive  Wiener  filter,  the  Widrow  filter  is  not  in  general 
optimum  even  asymptotically  [6]  ,  [71  .  However,  it  is  simple  and  in  a 
number  of  applications  it  performs  well  [8]  . 

2.  The  adaptive  frequency  domain  filter 

Our  objective  is  to  improve  the  design  of  adaptive  filters  in  the 
following  respects: 

1.  Reduction  of  the  number  of  adaptively  controlled  parameters. 

In  (1),  the  number  of  unknown  parameters  equals  the  length  N  of  the 
memory  of  the  filter.  This  constraint  complicates  in  many  cases  unnecessarily 


6 


the  design  of  the  filter.  We  mention  as  a  simple  illustration  the 
non-recursive  realization  of  a  one-pole  system.  In  this  case,  the 
coefficients  a^Cn]  approach  a  geometric  progression  involving 
a  single  parameter.  However,  if  zq  is  close  to  one,  the  required 
length  N  of  the  non-recursive  approximation  is  large. 

2.  Data  vector  with  nearly  diagonal  correlation  matrix. 

The  correlation  matrix  of  the  data  can  be  simply  inverted  if  it  is 
diagonal.  This  is  not,  in  general,  true  for  the  estimator  in  (1). 

3.  Design  based  on  the  local  spectral  properties  of  the  data. 

In  a  number  of  applications,  the  design  of  an  estimator  is  simplified 
if  its  adaptively  controlled  parameters  are  directly  related  to  the  frequency 
domain  characteristics  of  the  data.  For  example,  elimination  of  high 
frequency  noise  or  low  frequency  clutter  is  simply  achieved  in  the  frequency 
domain.  Again,  this  cannot  be  done  with  the  estimator  in  (1)  because  each 
frequency  component  of  the  output  depends  on  all  the  components  a^[n]  of  the 
vector  A[n]  .  Frequency  domain  processing,  however,  based  on  traditional 
methods,  involves  global  spectral  properties  and  is,  therefore,  unsuitable 
for  local  processing.  To  eliminate  this  difficulty,  we  shall  introduce  the 
concept  of  running  transforms. 

As  a  preparation,  we  impose  the  constraint  that  the  vector  A[nl  can 
be  written  as  a  product 

A[n]  =  WB[n]  (14) 

where  B  [n ]  is  a  vector  with  2M+1  components  b  [n]  and  W  is  a  2M+1 
by  N  matrix  to  be  determined. 
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Inserting  (14)  into  (2), we  obtain 

y[n]=  XT[n]  WB[n]  (15) 

With 

P[n]  =  WT  X  [n]  (16) 

(15)  yields 

y[n]  =  p"^[n]B[n]  (17) 

The  matTix  W  must  be  so  chosen  as  to  lead  to  a  small  value  of  M 
without  significant  increase  of  the  estimation  error  and  to  a  data  vector 
P[n]  whose  covariance  matrix  is  nearly  diagonal.  These  requirements 
can  be  satisfied  if  we  use  the  2M+1  significant  components  of  the 
Karhunen-Loeve  expansion  of  X  [n]  .  However,  the  underlying  transformation 
matrix  depends  on  the  statistics  of  X[n]  which  we  assumed  unknown. 
Furthermore,  in  this  approach  the  spectral  properties  of  x[n]  are  not 
utilized  directly. 

In  the  course  of  this  investigation,  it  will  be  shown  that  a  matrix  W 
whose  elements  w^^  are  the  roots  of  unity: 

w.  =  wkm  where  w  ^  ^  (18) 

km 

has  the  desired  properties. 
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The  elements  of  the  resulting  data  vector  P[n]  will  be  denoted 
by  X[n,  m]  .  Inserting  (IS)  into  (16),  we  obtain 

N-l 

X[n,m]  =  \  x[n-k]  w^m  (19) 

k=  0 

The  sequence  X[n, ml  so  formed  is  the  mth  discrete  Fourier 
series  [9]  (abbreviation:  DFS)  of  the  N  most  recent  values  of  the 
process  x[n]  .  It  describes  thus,  the  local  spectral  properties  of  x[n]. 

A 

The  output  y[n]  of  the  resulting  estimator  is  given  by  [see  (17)] 

M 

y[n]  =  ^  X[n,  m]bm[nl  (20) 

m=-M 

and  it  defines  the  proposed  frequency  domain  filter.  This  filter  consists 
of  a  memory- less  system  (Fig.l)  transforming  the  data  vector  X[n]  into 
the  DFS  vector  X[n,  m]  followed  by  2M+1  weights  b  [n]  and  an  adder. 

The  weights  b  [n]  can  be  adaptively  controlled  to  minimize  in  some  sense 
the  estimation  error. 

The  advantage  of  this  filter  is  the  fact  that  we  are  operating  directly 
on  the  local  frequency  components  of  the  input.  Furthermore,  as  we  shall 
show  [see  (41)  below]  ,  a  significant  data  reduction  results  because  in  many 
cases,  the  required  value  of  M  is  much  less  than  N.  Finally,  the  correlation 
matrix  of  X[n,m]  is  nearly  diagonal  [see  (42)  below]  . 

The  disadvantage  of  the  frequency  domain  filter  is  the  need  to  compute 
the  sum  in  (19).  This,  however,  is  not  a  serious  problem.  With  the 
available  FFT  algorithms,  the  computations  can  be  carried  out  routinely. 
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Furthermore,  it  is  possible  to  compute  X[n,  m]  recursively  using  one 
multiplication  only.  Indeed,  from  (19)  it  follows  readily  that 

X[n,  ml  =  wm  X[n-1,  ml  +  x[n]  -x[n-N]  (21) 

We  discuss  next  the  frequency  domain  form  of  the  estimators 
introduced  in  Sec.  1  . 

The  frequency  domain  Wiener  filter.  We  wish  to  determine  the  2M+1 

A 

components  bm[n]  of  the  vector  B[n]  such  that  if  the  sum  y[n]  in 

(20)  is  used  as  the  estimate  of  a  process  y[n]  ,  the  resulting  MS  error  is 

minimum. 

If  the  processes  x[n]  and  y[n]  are  jointly  stationary  with  known 
second  order  moments,  then  the  optimum  B[n]  is  independent  of  n  and 
it  is  given  by  (orthogonality  principle) 

B  =  R'1  r  (22) 

where  R  is  the  correlation  matrix  of  the  running  DFS  X[n,  m]  of  the 
data  x[n]  ,  and  r  Is  the  cross-correlation  vector  between  X[n,  m]  and  the 
signal  y[n]  to  be  estimated.  Thus,  the  elements  of  R  and  r  are  independent 
of  n  and  are  given  by 

E(X[n,m]X#[n,k]}  (23) 
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and 


E  [  X[n,  m]  y[n]}  (24) 

respectively. 

In  the  next  section,  we  express  R  and  r  in  terms  of  the  moments 
of  x[n]  and  y[n]  . 

We  note  that,  if  2M+1  =N  ,  then  the  estimator  (20)  obtained  with 
the  vector  B  in  (22)  is  identical  with  the  estimator  (1)  obtained  with  the 
vector  A  in  (6).  As  M  decreases,  the  design  is  simplified  without 
appreciable  increase  in  the  resulting  M.  S.  error.  In  fact,  if  the  process 
to  be  estimated  is  narrow-band,  then  a  nearly  optimum  estimator  requires 
a  small  number  of  terms  in  (20)  [see  (41)  below]  .  No  such  simplification  is 
possible  with  the  estimator  in  (1). 

We  now  assume  that  no  prior  statistics  are  known.  In  this  case,  the 
estimator  is  a  time- varying  system  and  the  vector  B[n]  is  given  by 

B[n]  =  R_1[n]  f[n]  (25) 

A 

In  the  above,  R[n]  is  a  matrix  whose  elements  rmfc[n]  are  computed 
from  the  recursion  equation  [see  (11)] 

rmk[“+I]  =  <*rmk[n]  +*( (1-a)  X[n,  m]  X*gi,  k]  (26) 

A 
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and  r[n]  is  a  vector  whose  elements  y  [n]  are  computed  from  the 
recursion  equation 


Ym C n+l]  =  aYm[n]  +  (1-3)  X[n,  m]  y[n]  ( 

Reasoning  as  in  (12),  we  can  show  that  the  vector  B[n]  satisfies 
the  recursion  equation 

B[n+1]  =  B[n]  t (1-CL )  R-1[n+l]  e[n]  D[n]  (; 

A 

where  e[n]  is  the  estimation  error  and  D[nl  is  a  vector  with  components 
X[n,  m]  . 

The  frequency  domain  Widrow  filter.  In  this  filter,  the  components  b  [n] 

-  m  L  J 

of  B[n]  are  such  that 

t>m[n+l]  =  bm[n]  +  iJme[n]  X[n,  m]  (2 

as  in  (13).  Unlike  (13),  the  coefficient  (j  in  (29)  is  different  for  each  m. 
This  is  so  because  it  might  be  desirable  to  have  a  different  adaptation  time 
for  each  frequency  component.  We  mention  as  example,  the  use  of  this 
filter  for  speech  processing  [10]  . 

3.  The  running  DFS 

The  running  DFS  of  a  discrete  signal  f[n]  is  the  sum 


F[n,m]  f[n-k]  w'1 


12 


where  N  is  a  given  integer.  For  a  given  n,  F[n,m]  is  the  DFS  in  the 
variable  k  of  the  segment  f[n-k]  of  f[n]  shown  in  Fig.  2. 
Properties.  Suppose,  first,  that  f[n]  is  a  deterministic  signal  with 
z-transform 


F(z)  =  ^  f[n]  z  n 


In  this  case,  the  sequence  F[n,  m]  has  a  z-transform  in  the 
variable  n: 


F(z,m)  =  ^  F[n,  m]  z  n 


From  (30)  and  (32)  it  follows  that 


m  ,  ,N 


T<z.m>>)  z'k  F(z)wfa"  -  F(z) 

'  I  «  HI  I 


1  m  / 
1-w  /  z 


because  the  z-transform  of  f[n-kl  equals  z  F(z). 

This  shows  that  the  sequence  F[n,  m]  is  the  output  of  a  linear 
system  with  input  f[n]  and  system  function 


J(z,  m) 


i  m  / 

1-w  /  z 


(Fig.  3a)  because  w1N  =1. 
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We,  next,  assume  that  f[n]  is  a  stationary  random  process  with 
autocorrelation 


R[k]  =  E  {  f  [n+k]  f[n]  } 


(35) 


and  power  spectrum 

m 

S(tti)  =  y  R[k]e"lkujT  (36) 

k5-oo 

In  this  case,  F[n,m]  is  alsoarandom  process,  stationary  in  the 
variable  n  and  non- stationary  in  the  variable  m.  We  shall  determine  its 
autocorrelation 


RFF^k,m3  =  EtFLn+k,m]  F  [n,m]] 


(37) 


and  power  spectrum 


S  p,p,  (uu ,  tn )  —  ^  Rp  p  L  1  ® 
k=-» 


■  jkuuT 


(38) 


Since  F[n,  m]  is  the  output  of  the  filter  of  Fig.  3a  with  input  f[n]  , 
it  follows  that  [9] 


SFFh,  m)  =  S(uj)  |  J(ejujT,  m)  |2  =  S(uu)  ^ 


a,2  N>»T 

sin 


sin  (■ 


(39) 


"W‘ 
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The  inverse  transform  of  the  above  yields 
N-l 

RFF[k,m]=\  (N-|r|)  R[r  +  k]wrm  (40) 

r=-  N+l 


Equation  (39)  justifies  our  claim  that  frequency  domain  filtering 
leads  to  data  reduction.  Indeed,  applying  the  inversion  formula,  we 
obtain 


U 

E{  |  F[n,m]  f  }  =  ^  J  S(ou) 


sin 


2  NxT 


-a 


sin 


TTm\ 

TT) 


d!D 


where 


a  = 


n 

T 


(41) 


This  shows  that,  if  the  bandwidth  of  S(uj)  equals  B,  then 

2 

E  [  |F[n,  m]  |  }  takes  significant  values  only  for  m  of  the  order  of  BN  la 
(Fig.  4).  If,  as  it  is  often  the  case,  f[n]  equals  the  samples  f(nT)  of  a 
continuous-time  process  f(t)  with  power  spectrum  Sc(uu)  and  T  is 
sufficiently  small,  then 


S(uu)  =  Sc  (uu )  for  |  u)  |  <  a  =  TT  / T 


(no  aliasing).  If,  therefore,  T  is  sufficiently  small,  then  the  number  of 
the  significant  components  of  F[n,  m]  is  small. 
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We  sh'll  next  show  that 


E  [  F[n,  mj]  F'  [n,  m2]] 


{m2-mj)/2  „  sinf%T  ,in  M 

= - t— -  l  S(uu) - — - dou  (42) 

'J  .  /uiT  nml\  .  /  u)T  ffm2\ 

-a  Sln  ("2 - ?7~  )  sin  ( ~2 - JT"  ) 

As  we  have  shown,  the  processes  F[n,mj]and  F[n,m2]  are  the  outputs 

of  the  two  systems  of  Fig.  3b  with  common  input  f[n]  and  frequency  responses 


respectively.  Therefore,  the  cross-power  spectrum  of  the  processes 
F[n,mj]  and  F[n,  m2]  equals  S(yu)  multiplied  by  the  product  of  the  first 
fraction  in  (43)  times  the  conjugate  of  the  second.  Equation  (42)  follows, 
thus,  from  the  inversion  formula. 

It  is  easy  to  see  that  the  fraction  in  the  integrand  of  (42)  decreases 
rapidly  as  |m^-m2j  increases.  This  justifies  our  second  claim  that  the 
processes  F[n,m^]  and  F[n,  m2]  are  loosely  correlated  if  ]m^-m2|  is 
large  compared  to  one. 

Filtering.  We  show  next  that  an  arbitrary  non-recursive  filter  of  length 
N  can  be  simply  realized  in  terms  of  the  running  DFS  F[n,  m]  of  the  input 
f[n]  .  As  a  preparation,  we  establish,  first,  the  inversion  formula  for 
running  transforms. 
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W  e  maintain  that 


N-l 

f[nl  =  4r  }  FCn'm]  (44) 

m=  0 

Proof.  Since  F[n,  m]  is  the  DFS  of  the  segment  f[n-k]  of  f[n]  ,  it 
follows  from  the  inversion  formula  for  DFS  that  [9] 

N-l 

f[n-k]  =  ^y  F[n,  m]w'km  O^k^N-l  (45) 

m  =0 

and  (44)  results  with  k=0. 

We  now  assume  that  f[n]  is  the  input  to  a  non-recursive  filter  with 
delta  response  h[n]  and  system  function 

N-l 

H(z)  =  £  h[n]  z"n  (46) 

n=  0 

We  shall  express  the  resulting  response 
N-l 

g[n]=  y  f[n-k]  h[k]  (47) 

k=  0 

in  terms  of  F[n,  m]  .  Inserting  (45)  into  (47)  and  changing  the  order  of 
summation,  we  obtain 

N-l  N-l 

g[n] =  TT  \  F[n,m]  ^  h[k]w_km 
m=  0  k=  0 
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Hence  [see  (46)] 

N-l 

g[n]=4r^  H(wm)  F[n,  m]  (43 

m=  0 

This  result  shows  that  an  arbitrary  filter  with  system  function  L(z) 
can  be  approximated  by  a  non-recursive  filter  H(z)  such  that 

H(eJU,T)  =  L(ejujT)  for  «,=  ^m/N 

It  is,  thus,  an  interpolating  filter  [12],  [13]  ,  interpolating  the  desired 
frequency  response  at  N  equidistant  points. 

The  running  low- pass  filter.  A  special  case  of  (44)  leads  to  a  simple 
design  of  a  low-pass,  high-pass,  or  band-pass  filter.  We  shall  discuss 
only  the  low-pass  case.  The  result  is  a  refinement  of  a  familiar  method 
of  filtering  (smoothing)  used  to  estimate  a  slowly  varying  signal  contaminated 
by  rapidly  varying  noise  [14]  . 

We  first  observe  that,  since 


H(wm  1N)  =  H(wm) 


F[n,  m+N]  =  F[n,m] 


the  limits  in  (48)  can  be  changed  appropriately. 

The  running  low-pass  filter  is  a  non-recursive  filter  such  that 


H(wm)  = 


jm|  <;M 


otherwise 
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We  shall  determine  its  delta  response  h[n]  and  frequency  response 
H(eJa;T).  Inserting  (49)  into  (48),  we  obtain 

M 


M 

V 

m=l 


g  [n]  =  ')  F[n,  m]  =  F[n,  01  +  -jq-  ^  Re  F[n,  ml 

m=-M 

From  the  above  and  (30)  it  follows  that 

N-l 


M  N-l 

i  r 


M 


g[n]  =  -fr)_  )  f[n-k]w  £[n-k]  y  w 

m=-M  k=0  k=0  m=  -M 


km 


(50) 


Hence, 


M 

»[“]  -  tt  l 

m=-M 


(M+l)n  -Mn 
nm  w  '  -w 


w 


N(l-wn) 


sin(M+-^-)  — ^ 

^T 


N  sin 


"rr 


(51) 


for  n  between  0  and  N-l  and  h[n]  =0  otherwise. 
The  system  function  is  given  by 


N-l 


M 


M 


H(z)  =  y  (tt  wimi)  z 


nm\  _-n  1 


n=  0  m=  -  M 


7T  \ 


1-z 


■N 


m=  -M 


,  m , 
1-w  /  z 


N 


because  w  =1. 

As  an  illustration  of  the  running  low-pass  filter,  we  mention  the 
adaptive  estimation  of  a  process  y[n]  in  terms  of  the  data 


(52) 


f  [n]  =  y[n]  +  v  [n] 


(53) 


containing  the  noise  process  v[n]  which  we  assume  white.  If  we  use  as  an 
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estimate  of  y[n]  the  output  g[n]  o£  the  low-pass  filter  with  M=0  ,  we 
obtain  [  see  (50)] 

N-l 

SCn]  =  -J-  F[n.  0]  =  3j-  y  f[n-k]  (54) 

k=0 

This  is  the  average  of  the  N  most  recent  values  of  the  data.  The 
resulting  MS  estimation  error  depends  on  N.  The  component  of  g[n]  due 
to  the  noise  decreases  if  N  is  large,  however,  the  distortion  of  y[n]  (bias) 
due  to  smoothing  increases.  To  minimize  the  MS  error,  we  must  vary  N 
adaptively  [14]  taking  into  consideration  the  local  properties  of  y[n]  and 
v[n]  .  The  running  low-pass  filter  permits  us  to  determine  these  properties 
and  to  very  the  effective  width  of  h[n]  by  varying  M.  Furthermore,  it 
improves  the  estimation  because  it  uses  not  only  the  average,  but  also  the 
first  M  harmonics  of  the  segment  f[n-k]  of  f[n]  . 
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7.4  Numerical  Results 


We  apply  the  adaptive  frequency  domain  filter  to  the  problem 
of  estimating  a  signal  y[n]  in  terms  of  the  data 

xGn]  =  y[n3  +  v(n] 

shown  in  Fig.  5a.  The  unknown  signal  yLn]  is  generated  as  the  output 
of  a  fifth  order  low-pass  digital  filter  (Butterworth)  driven  by  a  white 
noise  process  *[n]  (Fig.  5b).  The  noise  component  vfc]  is  white  and 
the  signal-to-noise  ratio  measured  as  the  ratio  of  the  corresponding 
energies  in  the  observation  interval  equals  2DB. 

A  _ 

The  output  ytn]  of  the  estimator  of  yin]  is  the  sum 
M 

y[n]=  7  X[n,m]bm[n] 
m=-M 

as  in  (20)  where  XLn,  m]  is  the  mth  running  DFS  obtained  with  FFT  size 
N=32.  The  value  of  M  is  determined  by  establishing  a  threshold  level 
T  and  accepting  only  the  components  of  XLn,  m]  that  exceed  T.  The 
value  of  T  is  chosen  either  by  prior  knowledge  of  the  r.  m.  s.  level 
of  the  noise  v[n]  or  by  a  measurement  of  this  level.  Since  viji]  is 
white  by  assumption,  it  suffices  to  measure  the  components  of  XLn,m] 
for  sufficiently  large  m.  The  results  of  the  measurement  are  shown  in 
Fig.  5. 

In  Fig.  7  we  show  the  value  of  M  so  obtained.  As  we  see  from  the 
figure,  after  some  initial  interval,  M  varies  between  1  and  4. 


(55) 
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The  coefficients  bmtn]  are  adaptively  controlled  as  in  (29)  and 
the  output  y[nl  of  the  filter  is  shown  in  Fig.  8. 

For  comparison,  we  show  in  Fig.  9  the  estimate  obtained  with  the 
time-domain  widrow  filter  (13).  As  we  see  from  the  figure,  the  estimate 
is  not  as  accurate  as  the  estimate  of  Fig.  8  although  the  number  of 
adaptively  controlled  coefficients  is  considerably  larger. 


Appendix  A 

We  wish  to  estimate  the  mean  of  a  discrete  random  process  x[n] 
in  terms  of  the  weighted  average 

n-1 

x[n]  =  (1-a)  S  a*  xpn-k]  0<a<l 

k=  0 

of  its  past  values  x[n-k]  . 

We  assume,  first,  that  the  process  x[n]  is  stationary  with  mean 
^  and  autocorrelation  R(k]  : 

E  {  x[n]  }  =  n  E  [  x[n+k]  x[n]  }  =  R(k] 

From  (A-l)  it  follows  readily  that 
E{;[n]]=  (l-an)n 

-n)2]  >  Til  y  cJk|  [l-a  ' ,k|)]  R  M 

k=-n+l 

Equations  (A-3)  and  (A-4)  show  that,  as  n  ,  the  mean  of  x[n] 
tends  to  the  mean  of  x[n]  and  the  variance  of  x[n]  tends  to 


k=-» 


Therefore,  if  R[k]  ir.  absolutely  summable,  i.  e. ,  if 


(A-l) 

(A-2) 

(A-3) 

(A-4) 

(A- 5) 
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I=£  |Rpc]|  <-  (A-6) 

k=  -  m 

then  the  variance  of  x[n]  tends  to  zero  as  n  -»  and  a-*l  . 

From  the  above  it  follows  that,  if  the  constant  a  and  the  integer 

n  are  such  that 
o 

2  n 

(1-Ot)  I  < <  2-n  and  a  °  <<1  (A-7) 

then  the  mean  of  x[n]  is  nearly  r|  and  its  variance  is  small  compared  to 
rf  .  Hence,  x[n]  is  a  satisfactory  estimator  of  the  mean  of  x[n]  provided 
that  n  is  larger  than  nQ  . 

We  note  that  the  time-average  x[n]  can  be  determined  recursively: 

x[n+l]  =  a  x[n]  +  (1-a)  x[n]  (A- 8) 

with  the  initial  condition  x[l]  =  x[l]  .  This  follows  readily  from  (A-l). 

We  conclude  with  the  observation  that  x[n]  can  be  used  to  estimate 
the  mean 


E(  x[n]}  =  r|[n] 

of  a  non- stationary  process  provided  that  r|[n]  is  nearly  constant  in  any 
interval  of  length  n  , 
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Figure  5(b) 
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